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ABSTRACT 


Three thermal recovery processes are discussed: 
(1) cyclic or intermittent steam injection, (2) steam 


displacement, and (3) hot-water injection. 


An approximate mathematical model describing 
a steam displacement process for a radial system was 
formulated but not solved. The system of five non- 
linear partial differential equations with variable 
coefficients describes the simultaneous flow of three 
fluids, oil, water, and steam under the influence of 
a pressure and thermal gradient. The heat losses to 
the overburden and underburden were included in the 


formulation. 


A second mathematical model was developed 
to describe a linear hot-water injection process. The 
simultaneous flow of two fluids, oil and water, under 
the influence of a pressure and a thermal gradient were 
considered. The liquids were treated as compressible 
and the heat losses to the overburden and underburden 


were neglected. 
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The hot-water injection model was solved using 
numerical methods. The alternating direction explicit 
procedure (ADEP) and the Newton-Raphson iterative 
procedure were used to solve the hydrodynamical equations. 
Crank-Nicolson finite difference equations and the square 
root method for inverting a matrix, extended to non- 
symmetric matrices, were employed in the solution of 


the thermal energy equation. 


Pressure, saturation and temperature distributions 
along with cumulative oil production histories and a 
material balance check on the results are presented. A 
discussion of convergence problems encountered with 


the computer program is included. 


The alternating direction explicit procedure 
(ADEP) yielded a stable, convergent, computationally 
and economically feasible solution to a one-dimensional, 


compressible, isothermal, two phase flow problem. 


Secondly, mathematical modelling of a 
thermally stimulated reservoir model is possible by the 


method outlined in this work and the results show that 
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increased recoveries and lower injection pressures 
are characteristic of highly temperature sensitive 
Oil viscosity models when subjected to thermal 


stimulation. 


Although a verification of the model was 
not possible the general shape of the pressure 
distribution, saturation distribution, temperature 
distribution and cumulative oil production history 
curves are representative of those presented in the 


literature. 
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INTRODUCTION 


The future demands for oil appear favorable. The 
economic law of supply and demand requires that an increased 
demand should be followed by an increased supply to 
stabilize the unfavorable balance. The problem of in- 
creasing the supply may be met in two ways by the petroleum 
industry; 1) increase exploration, and/or 2) improve 


recovery techniques of the existing oil reserves. 


As a result of alternative number 2), cold- 
water flooding became a widely accepted secondary recovery 
mechanism. There still remained, however, a class of oil 
reservoirs which would not readily yield to production by 
primary or secondary recovery techniques -- the heavy oil 
(high viscosity oil) reservoirs. Thermal recovery processes 
were thus initiated as a means for production of heavy oils 


with extremely temperature sensitive viscosities. 


One such thermal recovery process, steam 
injection, has aroused great interest and controversy. 
Successes and failures of pilot projects involving steam 
stimulation of an oil reservoir have been witnessed. The 
reasons for success or failure in most cases have not 


been published but it is believed that economics plays 
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a major role.in most of the failures. 


The development of a mathematical model to 
Simulate an oil reservoir undergoing a thermal change 
appears desirable. A brief review of available analytical 
solutions (!-15)for describing some aspects of the thermal 
recovery process for hot-water, hot-gas, and steam in- 
jection was presented by Flock et al (16), They point 
out that the analytical solutions generally describe either 
the thermal phenomena or the fluid flow aspect but not 
both. Davidson (17) states that the thermal and hydro 
dynamical phenomena should be analyzed simultaneously 
and that a rigorous treatment of heat losses is essential 


to developing a realistic model. 


To describe completely and accurately a 
reservoir system subjected to thermal stimulation, a 
complex mathematical model is required. The resulting 
partial differential equations which satisfy the model 
and treat the thermal and hydrodynamical phenomena and 
the heat: losses. in.a.rigorous. manner are highly non- 


linear. 
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Realizing the shortcomings of attempting to 
solve this complex problem by analytical means, recent 
authors (1 6717/18719,20,21/22) have adopted numerical 


techniques which are solvable on high-speed computers. 


This work was done with the following purposes 
in mind: 
1) to review the theory and literature on steam 
injection (three-phase flow) and hot-water injection 


(two-phase flow). 


2) to develop a mathematical model which would 
describe the hydrodynamic and thermal behavior of 
a reservoir system, and to discuss the validity of 
the major assumptions incorporated in the 


formulation of the model. 


3) to discuss various numerical techniques 
available to solve the aforementioned 
mathematical model, and thence to employ a 
method which will make the solution to the 
problem both economically and computationally 


feasible. 
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THEORY AND LITERATURE REVIEW -- STEAM INJECTION 


THREE PHASE FLOW 


The intermittent steam injection process has 
been studied by a number of authors (!6,17,18,21,22,23,24) , 
The physical interpretation of this process would be as 
follows. Steam is injected into the entire reservoir 
thickness or into a highly permeable fracture for a 
specified time (one to three weeks CT Upon completion 
of the injection period the well is shut in and the 
pressure begins to decline. This soaking period (two 
days to one week(!7)), allows time for the transfer of 
the latent heat of vaporization of the steam to the 
reservoir system. Having allowed sufficient time for 
the steam to condense, the well is open for production. 
Production continues until the producing rate which 
prevailed prior to steam stimulation is approached, or 
until the economics of production becomes unfavorable. 
At this time steam injection is. resumed and the cycle 
of injection, soaking and production is repeated; 


hence the name intermittent or cyclic steam injection. 


A complete analysis of the cyclic steam 


injection process therefore involves three stages: 


1) a description of the injection process, 
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2) a description of the soaking (condensation) 
period, and 


3) a back-production analysis. 


Davidson (17) developed a general model to 
describe the intermittent injection of steam. Realizing 
the uneconomic possibilities of attempting to solve the 
general model, he reduced the model to three one- 
dimensional simplified cases. A portion of this work 
was presented by Davidson et al (24). Flock et al (16) 
presented a two-dimensional mathematical model to describe 
the injection period and a second model for describing the 


back-flow period. The soaking period was not described. 


The need for an approximate but computationally 
feasible model inspired other workers (18721722) to 
develop models describing the "huff-and-puff", "push- 


BilwesOon ovCclroesteam=inj ection process, 


The major practical drawback to this thermal 
recovery mechanism is the loss of production during the 
injection and soaking periods. A steam displacement 


process whereby injection and production of fluids occur 
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in different wells rather than in the same well as is the 
case for cyclic steam injection would eliminate the 
necessity for lost production time which is inherent in 
the cyclic process. The models developed for the 
injection period of the cyclic process should describe 


equally well the steam displacement process. 


Figure 1 is a representation of the physical 
system. The mathematical model to be developed departs 
from the true physical problem due to the underlying 
assumptions incorporated. The limitations of the model 
will therefore depend on the seriousness of these 
assumptions and on any restrictions imposed on the 
solution to the problem due to the use of numerical 


techniques. 


The major assumptions incorporated in the 


steam injection mathematical model are listed. 
Liye A mathematical procedure may be used to 
describe the physical system. 
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DIAGRAM OF THE PHYSICAL SYSTEM--STEAM DISPLACEMENT 
FIGURE 1 


permeability. 
The reservoir thickness is a constant. 


All fluids exhibit newtonian behavior in the 
reservoir. That is, a linear relationship 
exists between the shearing stress and the 
velocity gradient, with viscosity as the 


proportionality constant. 


The fluid velocity may be represented by 
Darcy's Law and may be applied to each phase 


independently. 


Arising from the use of Darcy's Law, fluid 


interaction is not considered. 


Darcy's Law is valid for nonisothermal flow 
(see Appendix D). 


The relative permeability concept is valia (25), 


Capillary pressure effects are negligible. 
This is equivalent to saying the pressure 
in each phase is the same and is equal to 


the system pressure. 
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10. There is no free gas or gas in solution 


in the reservoir. 


Tine Oil production increases are due to viscosity 
reduction, thermal expansion and fluid 
displacement. This assumes that there are 
no volatile components in the reservoir 
fluid and therefore steam distillation is 


not considered as a driving mechanism. 
Ls The wellbore heat losses are negligible. 


#35 The use of an effective thermal conductivity 
in the energy equation is acceptable (see 


Appendix D). 


14, The reservoir rock is considered incompressible. 


The problem resolves itself into describing the 
Simultaneous flow of three fluids, oil, water, and steam, 
under the influence of a pressure gradient due to the 
injection of steam and a thermal gradient as a result of 
the steam injection temperature. A detailed derivation 


of the steam injection model is presented in Appendix A. 
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The continuity equation for a compressible fluid 


may be written as follows: 


Rate of change Net mass flux by |Rate of mass in- 
of mass = bulk transport + |crease due to 


per unit volume per unit volume vaporization per 
unit volume 


Net mass flux by 
molecular 

diffusion per 
unit volume 


or in mathematical notation, neglecting mass flux by mole- 


cular diffusion, the statement becomes 


Sg teen’ (6.5,) = Vieip. Wi + NG (T,P) (1) 


where the subscript (i) represents a particular phase. 
(see nomenclature for symbol definitions) 
Darcy's Law neglecting gravity effects (horizontal flow) 


and capillary forces may be written for each phase. 
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For a radial system, the divergence (V-> K) in cylandrireéal 
coordinates for horizontal, one-dimensional flow is 


represented by (26) 


Combining equations (3) and (2) with equation (1), the 


hydrodynamical equations are obtained for each phase. 


Oil Phase 
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From material balance considerations, 
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The thermal energy equation as stated by Bird, Stewart and 


Lightfoot (27)is, 


Rate of gain of Rate of internal Reversible rate of 
internal energy = energy input by + internal energy in- 
per unit volume conduction per crease per unit 
unit volume volume by compression 
irreversible rate of net gain due to 
‘i internal energy in- condensation per 
crease per unit volume unit volume 


by viscous dissipation 


or in mathematical terms, the statement may be written as, 


o—-=-V-. qi, - p(Vev)-t:¥ v + V(TOrP) (8) 


The internal energy U is not a convenient property for 
working petroleum engineering problems. Manipulation of 
the basic thermodynamic functions enables the transfor- 
mation of equation (8) to an equation involving the system 
temperature and heat capacity rather than the internal 


energy. 
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The heat losses to the overburden and underburden may be 


described mathematically in one dimension by Fouriers Law. 
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The theoretical model thus consists of five non-linear 
partial differential equations with variable coefficients. 
The independent variables are r, t, and z while the 


dependent variables are P, Tor T1> Sor Su and Soto A 


realistic set of boundary and initial conditions for a 


) 


constant mass injection rate of water (G) and steam (Gop 


would be: 


iiictaL Conditions 
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To complete the model the fluid densities must be presented 
in appropriate forms as functions of pressure and temper- 


ature. 


Equations (4) through (10) with appropriate 
boundary and initial conditions and equations of state 
comprise the mathematical model to describe a steam 
displacement process. A generalized mathematical model 
describing the thermal recovery of oil in linear systems 
was presented by Gottfried (*°), Six coupled, non-linear, 
partial differential equations with variable coefficients 
were required to describe the three-phase fluid flow model 
which included conduction-convection heat transfer, 
reaction between oxygen and oil and aqueous phase change. 
The analysis was applied to the detailed simulation of a 
forward-combustion laboratory tube experiment. Two to 
three hours of computer time were required to simulate 


twenty-four hours of real time. 


The similarity of Gottfried's model and the one 
presented in this work would suggest that similar computer 
time might be required to solve the problem and this was 


deemed uneconomical. 
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THEORY AND LITERATURE REVIEW - HOT-WATER INJECTION 


TWO PHASE FLOW 


Although the complex mathematical steam dis- 
placement model developed in the previous section and 
described fully in Appendix A will not be solved in this 
work, it will serve as a useful basis for simpler models 
of hot-water injection or hot, noncondensible gas 


injection. 


In treating hot-water injection or hot, 
noncondensible gas injection it is also desirable to 
have an accurate description of the thermal and hydro- 
dynamical phenomena occurring. The analytical solutions 
reviewed by Flock et al (16)as mentioned earlier, did 


not satisfy both requirements. 


Fournier (13) proposed a model for prediction of 
oil recovery by hot-water injection incorporating the 
effects of viscosity ratio reduction, thermal expansion 
and heat losses to the overburden and underburden for 


a radial system. He treats the fluids as incompressible. 


A representation of the physical system for the 
present work is shown in Figure 2. Figure 3(b) shows 


that the system is an "open-flow" boundary problem. 
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DIAGRAM OF THE PHYSICAL SYSTEM -- HOT WATER INJECTION 


PLGURES 2 
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DIAGRAM OF "NO-FLOW" AND "OPEN-FLOW" BOUNDARIES FOR RESERVOIR 
MODELS 


FIGURE 3 
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Figure 3(a) shows another possible model, the "no-flow" 
boundary problem. The mathematical development of the 
two problems differs but it will be interesting to note 
that the finite difference approximations to the partial 
differential equations are identical as shown in Appendix 


D. 


This study will deal with the "open-flow" 
boundary problem as illustrated in Figure 2 and Figure 
3(b) and therefore only its derivation will be considered 
in this section. Appendix D shows the derivation of the 
"no-flow" boundary problem and a comparison between the 
partial differential equations and the finite difference 


formulations for both models is illustrated. 


The major assumptions employed in the "open- 
flow" boundary problem derivation are the same as those 
presented for the steam displacement model with the 
additional assumption that the heat losses to the 


overburden and underburden may be neglected. 


In considering hot-water injection or hot, 


noncondensible gas injection the problem involves a 
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description of the simultaneous flow of two fluids, 
reservoir oil and injected fluid, under the influence 


of a pressure and thermal gradient. 


For a linear hot-water injection process the 
combined continuity-Darcy flow equations (4) and (5) 


assume the form; (see Appendix B for the complete 


derivation) 

Oil Phase 
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The thermal energy equation (9) neglecting the 
viscous dissipation terms (assuming high velocity gradients 
do not exist in the reservoir) and the reversible inter- 
change of mechanical and thermal energy due to compression 


or expansion (important for gases only) becomes: 


With the appropriate substitutions the thermal energy 


equation assumes its final form; 
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The dependent variables are P, T, Sor and s. and the 


independent variables are x and t. 


Finally, specification of the initial and boundary 
conditions to be imposed on the mathematical model for a 
constant mass injection rate of water (G.) completes the 


problem. 
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Initial Conditions 


Os 
So (x, 9) Soi 


Patcx, 0) = P. 


Boundary Conditions 
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Equations (11) through (15) with the necessary boundary and 


initial conditions and appropriate relations for the fluid 


densities as functions of pressure mathematically describe 
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the hot-water injection model neglecting heat losses to 
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DISCUSSION OF POSSIBLE METHODS OF SOLUTION 


OF THE HOT-WATER INJECTION MODEL 


The solution of the hot-water injection problem 
will require a numerical technique. As stated previously, 
any restrictions imposed on the solution to the problem 
due to the use of numerical techniques will result ina 


limitation of-the model. 


Assuming a solution to the mathematical problem 
exists and will satisfy the physical problem, a method of 


solution must be prescribed. 


The simultaneous solution of the continuity- 
flow equations (equations (11) and (12)) with the thermal 
energy equation (15) would be desirable. An attempt at 
this solution was unsuccessful, however, no conclusions 
were drawn from the failure to obtain a valid solution. 
Rather than pursue the problem of trying to solve three 
non-linear partial differential equations with variable 
coefficients simultaneously, a second method was 


prescribed. 


Since the continuity-flow equations are both 


parabolic, their solution should be possible by some 
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existing numerical technique. Their solution would 
yield a pressure and saturation distribution over a 
given time interval. This would require that all 
temperature dependent properties would have to be 
calculated using a known value of temperature. The 
calculated pressure and saturation distributions could 
be used in the thermal energy equation to solve for a 
new temperature distribution at the end of the time 


step. 


Calculate Temperature Dependent 
Properties in Equations (11) 
and (12) using an initial speci- 


fied temperature distribution 


Solve Equations (11) and (12) 


Simultaneously to obtain a 
Saturation and Pressure 


Distribution(S.,S_, and P) 


Solve Equation (15) for the 
Temperature Distribution (T) 
Using the Calculated Values of 


Pressure and Saturation 
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Calculate Temperature Dependent 
Properties in Equations (11) and 
(12) using the Calculated 


Temperature Distribution 


Assuming this prescribed procedure will yield 
a solution to the problem, a suitable numerical technique 
must be chosen to solve the continuity-flow equations 


Simultaneously. 


The first attempts to obtain solutions of the 
two-phase flow problems by numerical techniques could be 
attributed to Muskat (28), Terwilleger et al(22)and West, 


Garvin and Sheldon (3°), 


Douglas, Blair, and Wagner (31) determined the 
behavior of a linear waterflood employing a modified 
Crank-Nicolson difference equation to approximate their 


second order, quasilinear, parabolic differential 
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equation. The difference equation was implicit and 


second order correct in both time and distance. 


Peaceman and Rachford (32) discuss two types of 


difference equations. 


bg explicit difference equations, which are 
simple to solve, but which require an 
uneconomically large number of time steps 
of limited size. The forward-difference 
explicit (FDE) procedure fits this 


classification. 


2) implicit difference equations, which do not 
limit the time step but which require at 
each time step the solution by iteration of 
large sets of simultaneous equations. The 
backward-difference implicit (BDI) procedure 


fits thils classification. 


The authors present an alternating-direction implicit 
procedure (ADIP) which requires. the line-by-line solution 
of small sets of simultaneous equations that can be 


solved by a direct, non-iterative method and is stable 


bus tinifvfani e@aw nolsaups ap 


+ eh 


.sonet2 tb bs oma tod nt 


to asaqyt ows ssuoeib(S*)baottions Bee, ’ 


st 
| 
t 


exe odd yBhoL BipS nad lant 


na SsiLupert fotde doa a) ble 


eyete emit to ssdmya spiel’ yllsoime open 


gontexsttib-biswiot ont: este betis ee 
eitds e312 sxudevoxd! {gaay - gtoet 
Sosa 


ton ob dotdw ,anoissups oonsx032 ip 320) 


ry ia y¢ 


+s guiupes oidw suc qete amit ntl 


hak. 
to noigerett ya noisuloe ods ils meal 


} rote 


- 29 =| 


for any size time steps. Iterative techniques are also 


discussed, 


Douglas, Peaceman and Rachford (33) employed 
ADIP to solve a two-dimensional, two-phase, incompressible 
flow model which included the effects of gravity and 
capillarity in the formulation. Blair and Peaceman (3+) 
extended the work of Douglas, Peaceman, and Rachford to 
include the effect of compressibility of the fluids. 
They successfully predicted downdip displacement of oil 
by gas in a laboratory model by numerical solution of 
the two-dimensional, two-phase, fluid flow equations 


using ADIP., 


Gottfried, Guilinger, and Snyder (35) present 
two numerical methods (explicit and implicit) for 
solving the equations for one-dimensional multiphase 
flow in porous media. Variable physical properties 
are included in the formulation although not used for 
the problems solved. Gravity and capillary forces 
were neglected. The explicit and implicit procedures 
were essentially of about the same order accuracy, with 


the implicit method requiring 15 to 20% more machine 
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time but allowing use of a larger time step for a given 


mesh spacing. 


Quon et al(36)discuss five numerical methods 
for handling parabolic differential equations - the 
forward-difference explicit (FDE), the backward-difference 
implicit (BDI), the Crank-Nicolson implicit (CNI), the 
alternating-direction implicit procedure (ADIP) and the 
alternating-direction explicit procedure (ADEP). Their 
main suggestion was that ADIP and ADEP computational 
models would describe equally well a physical or 
mathematical problem but that ADEP would offer greater 
computational efficiency. A second demonstration of 
the use of ADEP(37) shows that it offers an advantage 
over the FDE procedure from a stability viewpoint, 
thus permitting larger time steps; and it offers an 
advantage over ADIP from a computational point of 
view since only one non-linear algebraic equation had 
to be solved at a time rather than a set of equations, 


representing a row or column of grid points. 


Coats and Terhune (38)compared the accuracy 


and computing speed of ADEP and ADIP in numerical 
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solution of reservoir fluid flow problems by analyzing 
published example applications. They concluded that 
ADIP yields significantly greater accuracy but requires 


about 60 percent more computing time than ADEP. 


Carter (3%)offers further discussion to the 
Coates and Terhune paper. His experience with ADIP 
and ADEP indicates: (1) that ADEP extends the range 
of conditions for which explicit methods may be used 
in practice to obtain gas flow solutions; (2) ADIP is 
a method of more general utility than ADEP; and (3) 
the additional storage and computer time required by 
ADIP over ADEP should normally be of little concern 


with present day computers. 


This study is concerned with a one-dimensional 
two-phase flow problem and involves the simultaneous 
solution of two non-linear parabolic partial 
differential equations with variable coefficients. 

The alternating-direction explicit procedure (ADEP) 
has not been applied to a one-dimensional problem of 


this type as far as the author knows and hence the 


é¢ 


paisyisns vd emo! dotg state bi lt — 


*, 


gssas pebulonos ir sna al 


. Ta0A nad omid paisugaino oom ta 


ow 
Y ae 
-~ 
*.. ® 


ont on noreenoarb  tSsd3 xu sori esa 
q1G4 ‘ddiw sonsixeqxs Bit segs sitire 
opnst sit ebnotxe FAdA $d (Ly ieee 
beer sd yam eborson: dtoilqne dota baareet” 
sit GIGA (S) seanotstuto 2 wold eer. nisddo’ « 


(€) bas ;GacA neds yo Et ese cexaop | 


yd bexttupst emis tedugqmoD bins apres 


mgonon sitadirl #0 da ult sinton stu 


énookatlowin stit eoviloynt ded 
isisisq sitodiasd, sedanonoe w3 
ejnabotatees oldsitey aviw aac 
(4aqaA) ewenndd. uote oop 


to cian: isaoiene 


t 


samiehiae) 


solution to this problem was carried out by ADEP along 


with the Newton-Raphson iterative technique. 


Crank-Nicolson finite difference equations and 
the. square root method for inverting a matrix, extended 
to non-symmetric matrices, were employed in the 


solution of the thermal energy equation. 


Appendix C shows the detailed discretization 
and final formulation of the finite difference equations 


which were solved on an IBM 7040 digital computer. 
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DISCUSSION OF CONVERGENCE AND STABILITY 


OF THE COMPUTER PROGRAM 


In the application of numerical techniques 
convergence may be defined as the approach of the 
approximate numerical solution to the exact solution 
as the number of increments employed increases and in 
the limit approaches infinity as the time step size 
approaches zero. This criteria is sufficient if an 
exact solution is known, however, in many cases as in 
this work the exact solution is not known and some 
other criteria for convergence must be established. 
For the problem considered in this work convergence 
was assumed when the results obtained from using two 
different time step sizes (one time step size being 
one half the value of the other) compared within the 
"acceptable limits" of the problem. "Acceptable 
limits" are self defined. and could be termed as the 
tolerance allowed between the two sets of answers 
obtained when the computer program is executed for 


the two different time step sizes. 


The program, although stable for large time 
steps, was not satisfactorily convergent. Also, -for 


large time steps an oscillation in the pressure 
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parameter was witnessed for small values of real time. 
This oscillation dampened out with increasing values 


of real times. 


To obtain a satisfactory convergent solution 
it was necessary to employ very small time steps, in 
fact, so small that it became uneconomical in that one 
hour. of computer time was required to simulate two 


hours.of real time. 


In order to remedy the situation, that is 
to make the ratio of simulated real time to computer 
time as large as possible, the time step size must 
be increased or the number of grid points must be 
requced . 4This 1s <contrary-to the first definition. of 
convergence which stated that the grid size spacing 
should be decreased. It is thus evident that the 
grid size spacing may be increased at the expense of 
decreasing the accuracy of the solution. A decision 
must be made as to whether or not the accuracy of the 
results may be given some flexibility so that an 
economical solution may be obtained. The desired 


harmony between the accuracy of the results and the 
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size of the allowable time step will yield an economical 


and computationally feasible model. 


The critical time step concept for the Forward- 
Difference Explicit (FDE) method is not too useful for 
non-linear problems of this type as the fluid properties 
and hence the coefficients of the equations are not 
constant which suggests that the critical time step is 


BLsa no. a consiant. 


The work presented in reference (35) served as 
a useful guide in determining an approximate value for 
grid size spacing and time step size and for the 
general "“trend"’of the results to be expected for two- 


phase isothermal flow. 
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DISCUSSION OF RESULTS 


If an isothermal fluid injection problem is 
being considered the actual amount of oil displace- 
ment and hence oil recovery during, the injection 
process will depend upon the relative ease with which 
the two fluids (oil and water) can move through the 
porous media. The readiness with which a fluid flows 
increases with its saturation in the porous sand. 

During the initial stages of displacement the oil 
saturation in the model is high and the displacing 

water saturation is low and therefore from a relative 
permeability standpoint the oil will flow readily in 
comparison with the water. As displacement continues 
the permeability of the sand to oil decreases while 

the permeability of the sand to water increases (Figure 
Fel, Appendix F) until eventually large. volumes of water 
will effect removal of only a small amount of additional 


oi ie 


When a reservoir undergoes thermal stimulation 
from steam injection, it is believed that viscosity 
reduction, thermal expansion and steam distillation 


are the primary mechanisms responsible for additional 
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Oil recovery. For the hot-water injection model 
considered in this work only the viscosity reduction 
mechanism was considered and both the oil and water 
viscosities were treated as temperature sensitive. 

The function definitions for water and oil viscosities 


as given in Appendix F are: 
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where the value of (v7) determines the sensitivity of 
the oil viscosity to temperature changes as shown in 


Figure (F-3) in Appendix F. 


Thus, for the isothermal problems considered 
in this work oil production will depend on the 
displacement process only, while for the nonisothermal 
problems the oil production will depend on the 
displacement process and the viscosity reduction 


mechanism. 
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Since the reservoir oil has no ability to 
expel itself, the oil recovery will depend upon the 
availability of the displacing water and on its 
accumulation during the displacement process. The 
availability of the displacing water may be specified 
but to determine the accumulation of water and hence 
the efficiency of displacement it is necessary to 
know the pressure, saturation and temperature conditions 
existing in the model at various times and positions 
along the length of the model. The results presented 
in this work are the pressure, saturation and temperature 
distributions necessary to predict the oil recovery 
for an isothermal or nonisothermal water injection 
problem and are shown in the figures at the end of this 
section. These figures were plotted from the computer 
output presented in Appendix G. Although it was not 
possible to have a check on the results or the model, 
the curves are representative of classical or typical 


shapes found in the literature (35,49), 


Figures 4, 5, and 6 show the pressure 


distribution along the total length of the linear model 
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for indicated values of grid size spacing (Ax) and time 
step size (At) for one particular value of real time 

and for isothermal and nonisothermal injection conditions. 
Three values of the oil viscosity sensitivity factor (Y) 
were used; 0.56, 1.12, and 1.50 as shown respectively on 


Figures 4,5, and 6. 


Figure 7 is a plot of the saturation conditions 
in the model which correspond with the pressure 
distributions shown on Figures 4, 5, and 6. These four 
curves are presented to show the error. induced in the 
results if the grid size spacing is increased to allow 
use of larger time steps and less grid points so as 
to reduce the total computing time requirements. The 
curves represented by the smallest grid size spacing 
and time step size are used as reference curves from 


which deviations of the other curves may be compared. 


It-is. evident that the choice of grid size 
spacing (Ax = 0.300 feet) is not the best one from the 
point of view of accuracy but it was believed that 
the combined harmony of accuracy and total computing 


time required to solve the problem makes this choice 
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the most economical and computationally feasible with 
the existing computing facilities. Guided by this 
criteria, the computer program was run for the following 
conditions for both isothermal and nonisothermal 


injection, 


TABLE I 
v Ax (feet) At (hours) 
OuS6 0.300 0o01L5 
l1al2 0.300 0.020 
ha5o H2300 Ondis 


Gottfried et al(35)show a family of pressure 
curves for a two-phase flow problem restricted to 
isothermal and incompressible systems. A sketch of 
their pressure profiles is shown in Figure 23 of this 


work. The following points should be observed. 


(1) The linear portion of the curves represents the 
pressure distribution ahead of the "flood front". 
For compressible systems this part of the curve 


should not be linear. 
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The non-linear parts of the curves represent the 
pressure distribution behind the "flood front", 


that is, the two phase flow region. 


As real time increases, the "flood front" advances 
and the pressure drop through the model decreases 
in the region behind the "flood front" since the 


resistance to flow decreases. 


Figures 8, 13, and 18 in this work represent 


the pressure distributions along the model length for 


isothermal and nonisothermal injection conditions. The 


curves exhibit similar "trends" to the results of 


Gottfried et a1 (35), The following observations may be 


made. 


(1) 


Although the system (oil and water) was treated 
as compressible the curves still remained essent- 
ially linear ahead of the "flood front". This 

is attributable to the low Povacuton rate (low 
pressure) and the very slight compressibility of 
liquids at low pressures. For high injection 


pressure and for gas flow it is believed that 
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compressibility will be important and hence the 


mathematical model was made general. 


The nonisothermal pressure gradients are lower 
than the isothermal pressure gradients for the 
Same value of real time at large values of real 
time. The thermal stimulation supposedly de- 
creases the flow resistance and thus the 
pressure drop decreases in the thermally 


stimulated regions. 


The curves do not converge to one common linear 
pressure distribution line ahead of the "flood 
front" for small values of real time. The 
error is caused by the choice of a large grid 
Spacing and time step which was discussed 
earlier. For the nonisothermal system 
conduction ahead of the "flood front" stimulates 
the region as discussed in point (2) and the 
pressure drop should not be along the same 

line if conduction effects are important. 

The errors induced, however, appear to be 


tolerable. 
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Figures 9, 14, and 19 represent the saturation 
distributions and Figures 11, 16, and 21 show the 
cumulative oil production histories for isothermal and 
nonisothermal flow conditions corresponding with the 


pressure distributions shown in Figures 8, 13, and 18. 


TABLE Ls 
Time Oil Recovery Maximum 
(Hrs.) (Eraction.O£L L0Oi1l.an Place) .Prontal Advance "A 
(Feet) 
36 QO. 94 - (0.2254 )* Dae Sb2egorDs) 0.56 
36 OF {0s Soy ASS sy Ol eee is Teel 2 
36 O5194. (0,184) Die ene) 150 


* The bracketed values are for nonisothermal conditions 
and the unbracketed values are for isothermal conditions. 


It would seem reasonable that the greater the 
sensitivity of the.oil viscosity: to temperature. changes 
the faster the "flood front" would advance and the 
greater the oil recovery should be. Figure 9 illustrates 
a faster advance of the "flood front" for the nonisothermal 


system and Figure 11 illustrates. the greater recovery for 
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the nonisothermal system. To be more specific, looking 
at. Table It’ fora value’ of oil sensitivity factor 
(Y=0.56), the frontal advance for the nonisothermal 
system is 2.55 feet while the advance for the isothermal 
system is only 2.25 feet... Similarly the. nonisothermal 
system depletes 21.4 percent of the original oil in 
place as compared with only 19.4 percent for the: iso= 
thermal problem. Since the curves on Figure 19 and 21 
are for a higher sensitive oil viscosity (vY=1.50) they 
should show an even more pronounced frontal advance and 
greater recoveries. Observation of the curves shows 

the reverse situation occurring at large values of 

real” time.” Looking at’ Table’ It for’ a~ value’ of: ¥=1°.50 
the maximum frontal advance of 2.25 feet is the same 

for both systems while the recovery of the nonisothermal 
problem (18.4%) is less than the recovery of the 
isothermal problem (19.4%). This would seem to 
indicate that the use of the large grid spacing is not 
acceptable or that it induces larger and larger errors 
as the fluid properties become more and more non-linear. 
Figures 14 and 16 show respectively the saturation 


distribution and the cumulative oil production for the 
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case where the oil sensitivity factor (vY ) equals 1.12. 
Again undesirable results are obtained as seen from the 
curves and from Table II. The maximum frontal advance 
is the same for both isothermal and nonisothermal flow 
(1.95 feet) and the production under isothermal 
conditions (17.4%) is greater than the production 
under nonisothermal conditions (15.9%). For this case 
the combined effects of large grid spacing, large time 
step and highly sensitive oil viscosity produced the 
errors in the results. For the isothermal systems 

the results should be identical for all problems but 
looking at Table II it may be seen that for /Y=0.56 and 
Y=1.50 this condition is satisfied but for V=1.12 the 
results are not identical. This may be explained by 
inspecting Table tT. For v=0.56 and vY=1.50 the grid 
spacing and time step sizes (At) are the same, For 
Y=1.12 the grid spacing is the same but the time step 
is larger. The following deductions may be stated. 

For a small value of the oil sensitivity factor (v= 
0.56) the use of the large grid spacing (Ax=0.300 feet) 


and a time step size of At=0.015 hours is acceptable. 
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As the oil sensitivity factor increases the use of large 
grid spacings and large time steps yields undesirable 


results. 


It will be interesting to note that the isothermal 
systems have only pressure as the driving force while 
the nonisothermal systems have pressure and temperature 
driving forces. It was mentioned earlier that the 
pressure gradients for the nonisothermal systems were 
smaller than for the corresponding isothermal systems 
and that as time increased the differences became more 
meonguncea. in actual fact,, for the conditions: of-this 
problem the water will not breakthrough the oulet end 
of the system for the nonisothermal systems. Figure 
8 may be used to provide an explanation. It is 
observable that the nonisothermal pressure gradients 
are approaching the initial pressure faster than the 
isothermal gradients. Eventually the pressure 
distribution will equal the initial pressure distri- 
bution and there will be no pressure driving force 
and no flow. It happens that for this problem the 


pressures equalized before breakthrough and the only 
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solution would be to increase the injection rate and 
thus the injection pressure. The computer output in 
Appendix G (Computer Output I) for the nonisothermal 
system and for a real time of forty-five hours shows 
the pressure equalization. Since this phenomena was 
not anticipated there was no lower bound placed on 
the pressure and the results show the inlet injection 
pressures being lower than the outlet pressure and 
hence reverse flow occurs. A lower bound on the 
pressure such that it cannot go below the value at 
the outlet end of the model should have been in- 


corporated in the program. 


Finally, Figures 12, 17, and 22 show the 
cumulative material balances for each flowing phase. 
The material balance should equal 1.0 and deviations 
from this value indicate possible errors in the 
method employed or in the choice of the grid spacing 


and time step. The curves are self explanatory. 


In summary, the alternating direction explicit 


procedure (ADEP), written for one-dimension, was 
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successfully applied to solve an isothermal, compressible, 
two phase flow problem. The thermal aspect of the 
problem was treated separately and the heat losses 

to the overburden and underburden were neglected. The 
author is the first to admit that the developed model 

is not highly sophisticated from the point of view of 
Simulating an actual thermal model but the experience 
and insight gained from this study leads the author 

to believe that the combined thermal and hydrodynamical 
aspect as well as a rigorous treatment of heat losses 

to the overburden and underburden could be considered 
Simultaneously employing the same approach used to solve 


the simple isothermal two phase flow problem. 


Although the method appeared to be 
unconditionally stable, convergence problems were 
encountered. The choice of grid size spacing and 
ae step size for non-linear problems of this type 
must be made on a trial-and-error basis or from 
previous experience (which the author cannot claim). 
The choice made in this work was such that an 


economical and computationally feasible model could be 
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developed. The economics of required computing time 
should not impose any restrictions for future workers 
as the IBM-7040 computing facilities at the University 
of Alberta are presently being replaced by an IBM- 
360-67 computer which is approximately six times 
faster. Smaller grid size spacing and smaller time 
step sizes may be used making the results more 
accurate and allowing calculation of the highly non- 


linear problems possible. 
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CONCLUSIONS 


Based on the work carried out in this study 


and the results obtained from specific example problems 


the following conclusions are stated: 


(1) 


(2) 


(3) 


(4) 


(5) 


Mathematical modelling of a thermally stimulated 
reservoir model is possible by the method outlined 


in this work. 


Thermal stimulation of reservoirs containing highly 
temperature sensitive oil viscosities increases the 
oil recovery and allows use of lower injection 
pressures as compared to the isothermal injection 


problem. 


Although verification of the model and the results 
was not possible, the results are representative 


of those found in the literature to date. 


The alternating direction explicit procedure may 
be used to calculate a one-dimensional, compressible, 


isothermal, two phase flow problem. 


The greater the non-linearity of the fluid 
properties the more restrictive the choice of 


grid size spacing and time step size becomes. 
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The use of small grid size spacings and time 
step sizes will tend to eliminate convergence 
problems and the increased computing time re- 
quired should not be restrictive with the new 


IBM-360-67 computing facilities. 
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(1) 


(2) 


(3) 


(4) 


RECOMMENDATIONS 


The model developed in this work could be extended 


to multiphase flow systems. 


The combined thermal and hydrodynamical aspect as 
well as a rigorous treatment of heat losses to the 
overburden and underburden could be considered 
Simultaneously employing the same approach used 

to solve the simpler isothermal two phase flow 
problem in this work; that is, ADEP and the 


Newton-Raphson iterative technique. 


Due to the recent controversy about the use of 
ADEP it would be desirable to employ other 
numerical techniques to solve this problem and 
compare the results obtained by the various 


methods. 


In this work the same relative permeability 
eerrelations’ for Ko and oy were used for 

the isothermal and nonisothermal systems. This 
may not be acceptable and one reason that 
actual recoveries may be higher for a thermal 
model is that the Ko and Lae data might be 


more favorable for the thermal system, and 
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hence work in this area might be considered. 


(5) _.The mobility of. the. fluids hoe ue and Sb Fe 
may be such that if the mobility of the water 
changes faster than the mobility of the oil 
during thermal stimulation, recoveries may 
actually be less for. the nonisothermal system 
as compared to the isothermal system. A 
comparison could be possible by using this 
work for both favorable and unfavorable 


mobility ratios to see if the phenomena 


could actually be observed and predicted. 
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NOMENCLATURE 


Cc = fluid compressibility, psia7! 

CCl a= fieldszunitdconversion constant = 2.637(107-") £t2. 
BY a /J0G «lit. 

CO = cold oil zone 

ce = specific heat at constant pressure, BTU/lbm°R 

‘een = specific heat at constant volume, BTU/lbm°R 

CW = cold water zone 

e = exponential = 2.171828 

7g = acceleration of dravity = 32.17 ft./sec? 

G. = mass. injection. rate.of water, 1bm./hr.ft2 

Got = Mass injection. rate.of steam, lbm./hr.ft? 

h = permeable sand thickness, ft. 

HW = hot water zone 

k = effective permeability, md. 

K = absolute permeability 

ce = effective SHEA conductivity, BTU/hr. it. & 

L = length of model, ft. 

P = absolute pressure, psia. 

q = Darcy's Law velocity, ft./hr. 

Q = volumetric flow rate, ft?/hr. 

a= heatuilux,) BRU/hE.£ts 


r = r-direction coordinate, ft. 


‘“sieq 1 a Jdtenon 
$42 (*-OL) Tea. = _dossano9 nekEseWNDD 


‘eA, 


elon Se ce 
0 eee 


~~} 
PSs 
i 


aE ,sweestq snsden09 35 teod ol asi 


A? md. EXOT Camu Lov saasenoo 35 sted 9 tks 


| | essiti.s = “tata us 
Soesa\.dt FLI.S& = oie te stoi 
+3.4d\.mdl ,tedew #0 ett _motsostitt 


$s. ad\ etl MB27 2 to S352 suena 


2 eaondtoids base { som oq 
4 yh J 


ETS ait tor 


bm .yvit rn 


A° 3%. 2d\0Ta 


* | 


% 


rj 


external reservoir radius, ft. 
wellbore radius, ft. 

fluid saturation, fraction 
steam zone 

time, hrs. 


temperature. °F, 


absolute temperature of the reservoir sand, 


absolute temperature of the overburden and 
underburden, °R 

PLULE*veLlociley; tee7 ut. 

vaporization term 

x-direction coordinate, ft. 

y-direction coordinate, ft. 


z-direction coordinate, ft. 


GREEK LETTERS 


thermal diffusivity = k_/#Cp, ft?/hr. 
forward difference 

§@-direction coordinate 

fluid viscosity,’ cp. 
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le) = density, lbm./ft? 

* = shearing stress 

) = effective porosity, -fraction 
SUBSCRIPTS 

e = effective 

i = a Fig Rh eG le: ie 

INJ = injection 

5 = jth grid intersection in the x-direction 
fe) = oil 

ol = Piivial ol conditions 

OB = overburden 

he = rock and relative 

st = steam 

W = water 

wc = connate water 

wi = initial water conditions 

x = x-direction 

y = y-direction 


Z = z-direction 
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SUPERSCRIPTS 

(yy) = (ny th time step n= 0), 1,82, 
(n+l) = (nti)st.time step, ne=.0,.1, 
yh = (r)th recurrence parameter fr 
(y+1) = (r#1) recurrence parameter 
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APPENDIX A 


DERIVATION OF THE STEAM DISPLACEMENT MATHEMATICAL MODEL 
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DERIVATION OF THE STEAM DISPLACEMENT 


MATHEMATICAL MODEL 


Basis: Unit Reservoir Volume 


The simultaneous flow of three fluids, oil, water, 
and steam may be described mathematically in terms of the 


continuity equation for each phase. 


The continuity equation for a compressible fluid 


may be written as follows: 


Rate of change Net mass flux by Rate of mass 
of mass =/bulk transport +/increase due to 
per unit volume per unit volume vaporization per 


unit volume 


Net mass flux by 
+|/molecular diffusion 
per unit volume 


or in mathematical notation, neglecting mass flux by molecular 
diffusion the three phase flow written for each phase 


separately is represented as: 
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To (refers to reservoir temperature) 


The velocity of each phase may be represented by Darcy's 
Law neglecting gravity effects (horizontal flow) and 


Capillary forces, 
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Since a radial system is to be considered, the divergence 


(veA) in cylindrical coordinates may be written as, 


Re 132 1 JA, oA, 
VA=>-— (rA_) + ——— + (A-7) 
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Combining equations (A-4), (A-5), and (A-6) with equations 
(A-1), (A-2), and (A-3), and applying equation (A-7) for 
horizontal, one-dimensional flow the hydrodynamical aspects 


of the model are developed for each phase. 
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The Thermal Energy Equation 


The thermal energy equation as stated by Bird, 


Stewart, and Lightfoot (!) is; 


Rate of gain! of Rate of internal Reversible rate of 
internal energy| = |energy input by |+ |internal energy in- 
per unit volume conduction per crease per unit 

unit volume volume by compression 

irreversible rate net gain due to 

nC of internal energy condensation 
increase per unit per unit volume 
volume by viscous 
dissipation 


In mathematical terms the eguation is written: 


DU 
=> 
9—=- Ve dy ~ P(Vev)- 
Dt 
The internal energy (U) is not a convenient property for 
practical engineering problems. Manipulation of the basic 
thermodynamic functions enables the transformation of 


equation (A-12) to an equation involving the system 


temperature and heat capacity rather than internal energy. 
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THERMODYNAMIC RELATIONS 


Tao dU Base dS —«Pdy (A-13) 
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“Pdv GAre="-StdT = Pdv (A-15) 
vdP  aG s= -S dT + vaP (A-16) 

U = internal energy 

H = enthalpy 

A = total work Helmholtz Function 

G = Gibb's Free Energy 

S = entropy 


T,P,v = absolute temperature, absolute pressure and specific 


volume 
Assume S = f(v,T) 
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then dS = |— av + |— aT 
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It is desirable to have equation (A-17) in terms of 
pressure (P) and temperature (T) rather than entropy (S) 
and specific volume (v). An expression for the derivative 


[25 | in terms of P and T may be obtained from a Maxwell 
dt 


Relation. Since equations (A-13) through (A-16) involve 
point (intensive) properties (independent of path) they 
are exact differentials and therefore of the general 


form, 
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Applying Green's Theorem to equation (A-15), the desired 


Maxwell Relation is obtained. 
as rc 9P 
(33) = Ea (A-18) 
ai V 
Substitute equation (A-18) into equation (A-17). 
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as = ES dv. + =z aT (A-19) 
V 
Substitution of equation (A-19) into equation (A-13) yields 


an expression for the internal energy as a function of 


pressure, temperature, and heat capacity. 
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The substantial derivative of U, multiplied by the density 
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The specific volume (v) may be written as 


But the continuity equation states, 
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Substitute equation (A-21) into equation (A-20) to obtain: 
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Substitution of equation (A-22) into equation (A-12) 
renders the final form of the thermal energy equation, 
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The mechanical energy equation as taken from 


Bird, Stewart and Lightfoot (is, 


p — |— |= P(v-v) -v> (Pv) + p(¥eg)-(Ve[t+ev])+(t2VV)  (A=24) 


A comparison of the mechanical energy equation (A-24) with 
the thermal energy equation (A-12) reveals the terms 
P(V+¥) and (t:VV) are common to both equations and appear 
with the opposite signs in the two equations. Therefore 
these terms describe the interconversion of mechanical and 
thermal energy. The term P (V°v) can be positive or 
negative, depending upon whether the fluid is expanding 


or contracting; hence it represents a reversible mode of 
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interchange. The term (-t:Vv) is always positive and 
therefore represents an irreversible degradation of 
mechanical to thermal energy. These viscous dissipation 
terms may usually be neglected except for systems with 


large velocity gradients (1), 


Heat Losses to Overburden and Underburden 


To complete the mathemcatical model the heat 
losses to the overburden and underburden may be specified 


by Fourier's Law written in one-dimension. 


a°T) i aT) 
= — (A-25) 


2 
OZ op ot 


where T, refers to overburden or underburden temperatures. 


(1) "Transport Phenomena", Bird, R.B., Steward, W.E., and 
Digntroot, ©,N., John Wiley and Sons, Inc., Sixth 


Printing, 1965: 
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APPENDIX - B 


DERIVATION OF THE HOT-WATER INJECTION MODEL 
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DERIVATION OF THE HOT-WATER INJECTION MODEL 


Basis: unit reservoir volume 


For a linear hot-water injection process the 
combined continuity-Darcy flow equations (A-8) and (A-9) 


assume the form: 
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Equation(A-11) reduces to 
Sa.+S = 1.0 (B-3) 


The densities of oil and water as functions of pressure 
may be written as: 
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*refers to a reference value. 
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Mathematical Manipulations 
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Differentiate equation (B-4) with respect to pressure, 
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** Dranchuk, P.M. and Quon, D.: "Analaysis of the Darcy- 
ConginuLity Bouation”, Petrol...Soc. CIM, 17th 
Annual Meeting, Edmonton, May (1966) 
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Equation (B-7) therefore reduces to 
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Substitute equations (B-8) and (B-10) into equation (B-1). 
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Since the function definitions for the density of oil and 
water are identical the above analysis for the oil phase 
is valid for the water phase and the transformation of 


equation (B-2) may be written directly. 
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The Thermal Energy Equation 


The thermal energy equation (A-23) neglecting viscous 


dissipation terms and the reversible interchange of 
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mechanical and thermal energy due to compression or 
expansion becomes: 
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In cartesian coordinates the substantial. derivative is: 
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Applying equation (B-15) to equation (B-14) for one- 


dimensional, horizontal flow yields, 


oT oT a 
oC en = - pc. Vy Ai - Ve. ae (B-16) 


The velocity ee is given by Darcy's Law based on the 
effective flow area 
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Combining equations (B-17), (B-18), and (B-19) and assuming 

Peat. = TY = T (see Appendix D) and 
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results in the final form of the thermat energy equation. 
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The heat losses to the overburden and underburden are 


neglected in this model. 
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DISCRETIZATION OF PARTIAL DIFFERENTIAL EQUATIONS AND METHODS 


OF SOLUTION 
ALTERNATING DIRECTION EXPLICIT PROCEDURE (ADEP ) ==) ONE= 
DIMENSION 
From Appendix B, the continuity-flow equations 
(B-11) and (B-13) are: 
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Approximations to the spatial derivatives 


from the use of Taylor's Series, 
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Consider any point j on a grid spacing as shown in Figure 


(C-1). The following expansions may be written to third 


order correctness. 


Wy 


! 


QOHTIM aus. 2 


enoit BHUpP®S we 


bs 


situplti ai mwode 


Apsatado 


ee 


2HOUTANIOS 


ad. wv 


wv) 
ow 


. ie 
Mt 7; | 


y 


si 


; y 
ee 


Br 
ee 


ie 
Rat 
Rot 


al 


ay 
re ‘oi 

ALTUAAa TS a J 

Wo ; 

i "7 j Py “ ) 
orruaoa: #¢ eu: " 
4h Ge 8); wae AY 4 , 
\ a : 7S Pik a 
i i ash r ate i” 

ye 


furit goo 


4 


at 


ee 


oo} 
36 
§ 
whe f > oe | <a 
: ‘ ~ aie 
2 ; 
a? 


2evistsvirtob Isitsqe 


; 2ot192 


(A) $ 
me # CED 
ic 
e ~ 
25 


a, oR 


fuk 
. euttititjeaath 72 
| £99 

an 
| Q 
eee zs 
4, fOr 


@ Pach 20 98 


(x A 
sos 

‘ 
— 


$s 


(x) ruil uk” 


—_—— FORWARD SWEEP 


TIME (t) 
— REVERSE SWEEP 


=—o= 


ONE-DIMENSIONAL ALTERNATING DIRECTION EXPLICIT PROCEDURE -- 


GRID POINTS AND TIME LEVELS FOR THE FORWARD AND REVERSE SWEEP. 


FIGURE C-1 
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By adding equation (C-3) to equation (C-5) it can be 


shown that 
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By multiplying equation (C-5) by 8 and adding it to equation 
(C-4) an expression for the first derivative of T with 


respect to x is obtained. 
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equation (C-6) yields an expression for the first derivative 


(£5) My 


E(xaS) 
b=D) i 
( i 4 
(g=+9) “? er 
© (xas) 
(9=9) tg ——— 


i | 
1 
— 


~~ 


i | aval 7A 
noissups od tf eni bbe bts. 8 ya Daees > emi. 


oe \ : 

ee 

dtiw T to ov haa tem dexit oma ans ewes’ 
ae ; a NPA Si ke 


aR teic Oo gh ; +4 i $2" 
a aN ‘65 4 J . 


A ; ie = -, =v, 
a 4 } r ; 
’ rep. 1 A 
— bh 
eae Tae ie 4 ; 


/ 
i 


(8-5) 


hata seat 


ryt ‘ oa * A A) a A 
vw Nt nee mer 
ase joe 
i. | Pie wa ng 


-  105-- 


of T with respect to x incorporating two backward points 


and one forward point. 
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Equations (C-7) to (C-9) are the approximations to the 
spatial derivatives of pressure and temperature to be 


used in this study. 


Discretization of the hydrodynamic equations 
(C-1) and (C-2) using equation (C-7) yields equations 


(C-10) and (C-11). 
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At this point, the particular method used will dictate 
the time levels at which the various terms must be calculated. 
The alternating direction explicit procedure (ADEP) in one- 
dimension as shown schematically in Figure (C-1) was applied 
to equations (C-10) and (C-11). The boundaries at x = O and 
x = L are specified and it is desired to obtain a pressure, 


saturation and temperature distribution at all grid points. 


GRID POINT (j=2) - FORWARD SWEEP - FIGURE (C-1l-a) 
Oil 
k (n) 5p (n) b ce 5p (ni) 
oi a. ae oe 
u ox u ox 
ore es fe) 
TL 2 el 72 j-1 j-1 
Ax 
> 
3. ees g (ntl) 2 tae if NS Ute Ouse ig xP) (C-12) 
COMA Eile) : 4 J j 
k oP Go 
eee — = — = 0.0 (Boundary Condition) (C-13) 
U ox Cel. 4p 
= j-1 j-1 Ol x=0 
(n) 2 
k > (Ax) 
mc) (P..1-P.) wie ———[§0t8) si) ac sin?) (P De) oo=aa 
J J “(cot ay oe ae j j - - 


atetotb LLiw beam boddom: sspsnabiaiiee a 
betsivoiss ed seum saries ba i vauiaas 
-sro mi (@HGA) exubeso7g $i: 3 
beifags saw (l=) omupit nb ysteo: | 

bas 0 =x 3s esixsbavod oat tth-2) bad fers 7 
sisegstg & aistdo od benieeb et ti ‘bast ae 08 


,eatatog Bryp fis S65 


a \ mn Tr A Oa? 
wb mo 134 i je) L GH SS Woe 
— > 192 ie) ree, Se eee 


my PAY ge (EH gy (S\Eeet) tw 
{[~9) | ae | e 7) oe? +- | 


(El-0) (noftibaod ys 


= L097 = 


Water 
k (n) 5 P (n) k (n+1) Pp (n+l) 
Patisch an ei ee mi 
a aX ee OX 
jr L/2 ce aay2 j-1 4S 
AX 
) 
- —— | s{”)-s; 
CCL At aan) j 
k dP G 
Se pn = a (Boundary Condition) 
We ox ed Sail Bs 
jet aa) 
eRe (n+1) k (n) 
W i (n) 
CCl p 2 J 
- j+1/2 
> (Ax) @ 
4 ae Y g (ntl) 
COT (AC) 5 i 


+C (1-S 
W 


O 


(nl / a) 


j 


Oo. 


J 


es Webece(i-scmaayap en? 


(n+1) 
We 


meh era 
(C-16) 
Pf ean. 


(a) 4 
Bee =f a 


£-S) 


i gta se-) - 


| Dart) oo, , (S\ Dee) Lee 
weet Poe 


| ii | ‘ i | ra : 
5 re, a 
Ras 


“ 


C 


(qoitibaod yisbavod) 


é Ms noeve: io ‘ i ' - 
iP Sega: pen - ‘ 


ih pa ra 


P . : : ; , f rahe i > 
L om - fl i 7 : H i hits 


C 


H a t 1 aa 
a hg: J) CNet xs 
LR a Vhs 4 f/m 
A | ei © pe 
ee i ; : ° F 


d 
r 


- 108 - 


The following relations are applicable in all finite 


difference formulations in this work, 
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INTERIOR GRID POINTS (2<j< (L-Ax/2))-FORWARD SWEEP-FIGURE (C-1-b) 


Oil 
(n+1) (n) 
k k 
= a sala (P,P) ™ 
O Ho : : 
a et j+1/2 
o (Ax) 2 
lore eee g (ntl) a g (n) +C¢ g (nt1/2) (5 (n+1)_,(n) ) (C-25) 
GGM(at).| °3 °5 Sos J J 
Water 
(n+1) (n) 
k k 
WwW (n+1) WwW (n) 
- (P,_4-P.) = (Pea P.) 
Uy - Uw le 
jelLe 22 
_ $ (Ax)? i (n+1) (nt+1/2) (n+1) — (n) 
Ss = s poie( lees hae =P (C-26) 
ccl(at)| °3 °5 Wg (05 : te 


INTERIOR GRID POINTS (2<j< (L-Ax/2) )-REVERSE SWEEP-FIGURE (C-l-e) 


Ole 
(n) r (n+l) 
k k 
we) cll aie aa (n) 4. my (P.)-P.) (n+1) 
th J J Mo | J J 
j-1/2 ihe 


> (Ax) 2 
x pt mis (ee eth (2) ee ae ens (C-27) 
cel(at) | ° ty eta, ; ; 


(d-{-9) JAVOLT-IsIIWa | anAWsons tS 
cca nL ALL LLL LA 


(a) gg ¢) i im | 

t [+t ) bas F 
NEE \ a 

: | | ba ih. 


| Lid pie 
pg i ) (S\Ptohy, 5 + Me _ 
L 3 


(eS=P) i? oO ean ee 


ra 


lier Sails FE) 42+ ie 
4 re i 


(ae-9) aa 
, ye) 


/ ui 4 
io ie | s ig 


’ 
T ~~ | AO 
it 


Wi) roa ara? ee 


ee 
sap. 
SS sae 


ig \ 
a 


(9-[-D) SAUDI I-IaawWe desavasi= ((S\x4~7) of 2 gee 


a 
ive 
_ hd 


~ 


l i My, i" 
5 St 7 es ’ 
o> Pa : er , 


Ae lc ee yt Sy oP May ah Pee 
THR ae I ea Ae 
“roll a ettae * Rie) yy me? vl % 


= Lhb - 


Water 

(n) (n+1) 

k k 

W (n) W (n+1) 

Wy f ew 

j-1/2 j+1/2 
o (Ax) 2 | | 

a g (My o(mtlhis ¢ cg 2421/2) op (tL) _p (n) | (C-28) 

ccl(At)| °% ©; leg Bd J 4 


GRID POINT (j=L=-Ax/2)-FORWARD SWEEP-FIGURE (C-1l-c) 


Oil 
k (n) )P (n) k (n+1) Pp (n+1) 
0 indies ply § ead ger: 
U ox Wo ox 
Pat /2 41/2 je 1/2 4-1/2 
3/4 (AX) 
> 
a Tea th) a tia) at) Gy} neem 
ccl(at) | °3 °5 Pe 5 J J 
(n) (n+1) 
8 (k 4 (k 
pr (P..4-P ) (n) bh eae hig (P PD.) (n+1) 
8 thn 5 j 3 uy j 
> ERV/2 j-1/Y2 


> (Ax) 2 
., Sa sy? g (n+l) rr g (Nn) 6 alter Beate) 0) | (C-30) 


CCl (At) ©; °5 aes J 


({£+02) ti. , ) 
te mee is 


(Peet) ce A SRERAD el i 
ae ay ( O° ike 


sean | 


(o-{-D) SAUDE T-WaawsS Ca AWAOt= (S\x d= bat th 2 


——— 


(E46) 3 be sli . (a) 


4xe: 
O\ LE. §\lfec 


-_ -y - . ; ? tem 
(@$=-9) x 3 ) ie i is & ay. 
“Uea ey 

({+.1) 7 eee Se 


wet: ee 


; 6 
y UJ § 
ba 


at: 


if 


r - 
i+n),. s\n) ‘Ma, 
(O£~D) tae ie alii: aps ae 
| Fag neg eee 6 


= 122 = 


Water 
(n) , (n+1) 
8 {k 4-{k 
_ |_w PP een \giT) oan | 10, ioe he eS 
3= (hy jtl “3 cut j-l “3 
Po +1 /2 “3 4a /2 


9 (Ax) * (n) (n+1) ‘Geo, eran. 
gM) _ g + Cle (P| -p))| (¢-31) 
CC1 (At) j j °5 J | 


GRID POINT (j=L-Ax/2)-REVERSE SWEEP-FIGURE (C-1-d) 


Oil 
k a ) Dp (n+1) k (n) ) Pp (n) 
Ves - | ps0: px" 
u ox u 0X 
ea 1/2 j+1/2 ES Le) ee ae 
3/4 (Ax) 
a $ lay a THe temtenL/ 20 eFC eae 
CCl (At) 25 5 2 194 7 j 
(n+1) (n) 
8 fk 4 (k 
bel 22 @ to splywS the eaualtiShs eax iP.) 
3 O18 j+l 3 sult Tse 
0 441/72 Rirj- V2 


(Ax) 2 | 
est lane (Mhoy ag tg (PtAA2) pp MEbLerpA2d) (C~33) 
ccl(at) | °3 34 cee : ; 


F , een G ( 
(£-9) hag in ale t+) on) oy * 


C 


‘ y ay 
(pe f-0). SAUDTI~Saswe Bepavads(S\t ards EhS 


x6 ” ye ie 

cyst S\d-t 9 SNRRE | ae 
; ve a! nals “ : 
Cent ONE) ib ee 
; ai ae yi ' 
: J re ‘(ia 7 

bt) | (P+) a, (S\Eta) (oh, -. (tag 
a 380 + ~ a: 5, 


Ae: 0 le 


(SE-9) ( : es eg 


Shite! 


1h sags Re) ee anna al 
(ita), (S\ita), 4 4 (0) g thre. 
(€£=9) ims -{ eae nh i 


ees 
1 y , } 
’ eA en ( 
oa) apy 4 
| ara ps 
tt Oe a ae 
‘ / : vA 1,4 Via i iy oF we) i rae t i . Bi be ee, 
=f ! i ; AA wy { ats y aris } i tie ‘ ; ‘ { ¢ ; 
° eye Waa ts eae 4 (alee uM , 


F ‘ ie 0 is : } i 5 one Lie ' | 
ihcy Si liei l 
‘ Vi P al 2 7 ; Fj a y | 
‘ ese sedi it Ps - ' 


“ EES * 


Water 
/ (n+1) (n) 
8 fk Ciel ic 
ae acd (Pe =p ) (n+1) het ayy. (Pp. =pe) (n) 
ee RE j+l “3 nat 2 hh ae 
Oi /] 2 Lass eh iy do) 


(Ax) * 
ES ere MELD Oe ore PGUEUN IES (C-34) 


CC1 (At) 5 of ae j 


Newton-Raphson Iterative Procedure 


If the solution to the problem is started at the 


first grid point beyond the inlet boundary, the only unknown 


variables in the finite difference formulations are Se 
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and panic! and any fluid properties which are functions of 


these two parameters. Since the fluid properties are 
dependent on the two variables an iterative technique will 
be required to solve the equations simultaneously. The 
Newton-Raphson iterative procedure should suffice. Although 
the method does not guarantee convergence, especially for 
non-linear problems, no difficulty is expected as there 

are only two variables and an initial guess of these two 


parameters can be made reasonably close to the physical 
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The solution of equation (C-35) will be the next, supposedly 
improved, guess. 
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A recurrence relation for the (r+l) iteration may be 
formulated as follows; 
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where Es is specified-sufficiently: small to satisfy the 
problems's required accuracy. 

It 1s easily seen that the inversion of the 
Jacobian Matrix will dictate the ease of solving the 
problem by this technique. The problem, with only two 
Variables, will require the solution of a 2.x 2 matrix 
and this is very straightforward and should offer no 


Gerrveulty . 


To be more specific, if the iteration variables 


pei) tls) a Ate relhy) 
are denoted by Pe = PS and Sox = So. 


i) 
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It should be noted that the method is first order correct 
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The Thermal Energy Equation - Crank-Nicolson Approximation 


Equations 


(B-20) and 


(B-21) in Appendix B describe 


the thermal phenomena of conduction and convection in the 


model. 
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is a function of saturation only. 
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Using the Crank-Nicolson procedure, the finite 
difference approximation to equation (C-50) will be 
implicit and hence a series of equations must be developed 
across the grid network of Figure (C-1) until the outlet 
boundary is reached. Since the outlet boundary condition 
is specified the problem may now be solved. The solution 
involves the inversion of a matrix and it is assumed 


that this inverse exists. 


Equation (C-50) in discretized form using 
equations (C-7) and (C-8) to approximate the derivatives 
may be used to describe the inlet boundary point and 


interior grid points. 
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Forms of equations (C-7) and (C-9) to approximate the 
derivatives are required to describe the last grid point 
before the outlet boundary. 


The discretized form of the equation is: 
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Applying the Crank-Nicolson procedure to equation (C-52) 
results in: 
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Rearranging equation (C-54) so that terms to be evaluated 
at time (n+l) are on one side of the equation and terms 
to be evaluated at time (n) are on the other side of the 
equation, the final formulation of the finite difference 
approximation to equation (C-50) for the inlet boundary 
and interior grid points is developed. 
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Similarly, applying the Crank-Nicolson procedure to 


equation (C-53) results in: 
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Rearranging equation (C-56) in a like manner to equation 
(C-55) the final formulation of the finite difference 
approximation to equation (C-50) for the outlet end of 


the model is developed. 
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The use of third order correct approximations for 
the derivatives results in the formation of a non-symmetric 
matrix of equations. The square root method, extended to 


non-symmetric matrices, was employed to solve the problem, 
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APPENDIX D 


COMPARISON OF THE "OPEN-FLOW" AND "NO-FLOW" 


BOUNDARY PROBLEMS 
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Comparison of the "Open-Flow" and "No-Flow" Boundary 


Problems 


Figure 3 illustrated two linear models: (a) "no- 
flow" boundary case and, (b) "open-flow" boundary case. 
The primary difference in the mathematical treatment of 
the two models for an injection or production problem 
lies in the specification.of the boundary conditions. 

For model (a), there is no flow across the inlet and 
outlet boundaries. .A source term is required to describe 
the injection into or production from the model but 
without the necessity of specifying where the source 

term is located. Model (b) on the other hand specifies 
that flow may occur across the inlet and outlet boundaries 
and hence fluid may be injected or produced without 
specifying a source term. Model (b) was discussed in the 
section of "Theory and Literature Review - Hot-Water 
Injection - Two Phase Flow" and was developed fully in 


Appendix B. 


Model (a) will be developed in this section. 
The continuity-Darcy flow equations may be written for 
each phase with the mass injection rate of water being 


represented by a source term, 
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CONSTANT (at injection point) 


0 (at.all other grid points) 


Since the reservoir thickness, h, is not a function of x 
for this problem and the mass-rate input of oil is zero 


equations (D-1) and (D-2) may be written as: 
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Equations (D-3) and (D-4) may be manipulated in the same 
manner as the hydrodynamical equations (B-1) and (B-2) 


were in Appendix B giving rise to equations (D-5) and 


(D6) . 
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Equations (B-11) and (B-13) in Appendix B are the 
analogous hydrodynamical equations for model (b), the 


"open-flow" boundary case. 
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The source term will be equal to zero for interior grid 


elements and the two sets of equations become identical. 


Equations (D-5) and (D-6) will be discretized 
for the inlet boundary (j=2, Figure C-1). The approxi- 
mations for the derivatives will be the same as those 
used in Appendix C and the one-dimensional ADEP will be 


applied. 


Oil Phase (Forward Sweep) 
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C which shows that the fluid flow equation for the oil 


(8-G) 


ee (ft) g 


(@-a) 


" phbanags: ves) noise 


bixe torresad 10% oxes OF) taupe 


ae I He a 
.isoisaebi emooed en0i38upe to-etoe 7 be ae ‘ 
ie 


-ixouggs ofT - (49 owen set) ysbe wod ont ‘ 


J ne : ds 
pesitexosib ed iLtw. (are, ‘Das tend Boat 


aeons es 9ise ods ad Loew eavisavine 


~ 
a 
at 


ed ILiw WGA renvtbertoebieeans ant brs 9 


(Li) jee CLOT 
a6 aj 


— 


| yo \e6 ye 
46 528 ae a f° 
a 8.0 +i es 
36 6 


a de 


et ge SP oti 
a 


© 


{io ens aaa 


- 130 - 


phase in finite difference form is the same for both 


model (a) and model (b). 


Water Phase (Forward Sweep) 
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Equation (D-11) is identical to equation (C-17) 


in Appendix C which shows that the hydrodynamic equation 
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for the water phase in finite difference form is the 


same for both model (a) and model (b). 


In summary, although the mathematical development 
of the partial differential equations to describe the 
hydrodynamical aspect of the problem is different, as 
shown by a comparison of equations (D-5) and (D-6) for 
model (a) with equations (B-1l1) and (B-13) for model (b), 
the resulting finite difference formulations are identical 
as shown by a comparison of equations (D-19) and (D-11) 
for model (a) with equations (C-14) and (C-17) for model 
(b). This suggests that the solution to the problem 
will satisfy either the "no-flow" boundary case or 
the "open-flow" boundary case if the discretization of 
the partial differential equations is carried out as 


described. 
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(a) Discussion of Validity of Darcy's Law for Non-Isothermal 


Systems 


Darcy's Law for a newtonian liquid of constant density 


or of density as a function of pressure only may. be written 


as 
> 
gq = oo Vo 
ok 
where oO = — 
u 
Pp 
and “) = gz - — = POTENTIAL 


Since ® is a scalar potential and gq vsaadfluxkvector it is 
essential that gq be irrotational. This implies that the 


circulation (Ir) must equal zero, where 


$ - refers to a surface 
Hubbert (!) indicates that this equation will be valid for 
isothermal systems. Davidson?) shows for nonisothermal 
reservoirs the circulation is not zero and hence a rotation- 


al flow field develops. Strictly speaking it appears 
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Darcy's law is not valid for nonisothermal systems. 
However, the magnitude of the deviation in flow behavior 
if specified by Darcy's Law under nonisothermal conditions 
is unknown. Lacking an alternative, this work assumes the 
effect of nonisothermal conditions on Darcy's Law to be 


negligible. 


(b) Validating the Use of an Effective Thermal Conductivity 


in The Thermal Energy Equation 


The use of an effective thermal conductivity appears 
to be an acceptable approach in describing heat transfer 
problems using Fourier's law. Bailey and Larkin/(3) 
demonstrate the assumption of instantaneous heat transfer 
between the solid and liquid phases of the porous system 
is quite reasonable. This work therefore assumes the use 
of an effective thermal conductivity is valid and hence the 
temperature of the reservoir rock and its. surrounding fluid 


are the same, 
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APPENDIX. «F 


RESERVOIR AND FLUID PROPERTIES - COMPUTER INPUT PARAMETERS 
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RESERVOIR AND FLUID PROPERTIES - COMPUTER INPUT PARAMETERS 


) = Porosity = 0.24 

ep = Rock Density = 165.2 lbm/ft? 

ps = Oil Density at (P*=14.7 psia) = 55.0 lbm/ft 3 

px = Water Density at (P*=14.7 psia) = 62.4 lbm/ft ° 

us = Oil» Viscosity at 60°F = 265 cp. 

EPS. = Specific Heat of Rock = 0.23 BTU/1lbm°F 

Cp. = Specific Heat of Water = 1.0 BTU/lbm°F 

co. = Specific Heat of Oil = 0.650 BTU/lbm°F 

K = Absolute Permeability = 1000 md 

ke = Effective Thermal Conductivity = 0.5 BTU/hr.ft°F 
ea = Oil Compressibility = 3.0(10-°) psia@} 

ee Water,Compressib@lity =. 3.2({10%°) psiaq} 

Sic = Connate Water Saturation = 0.20 

Soa = Initial Water Saturation = 0.25 

Soi = pavciak Ort Saturation "=70..15 

P; = Initial Pressure Distribution = Constant = 14.7 psia 
T. = Initial Temperature Distribution = Constant = 60°F 
T = Water Injection Temperature = 200°F 
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FUNCTION DEFINITIONS 


=—_p*x 
Cc, (P px} 
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Oil Density BAR). Se 5 
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Water Density p, (P) 


Effective Permeability to Oil (Figure F-1) 
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Effective Permeability to Water (Figure F-1) 


Ere or Pic } 
k(S3) = K ma On 0 ie Se 0.8 (1) 
1-S 
wc 
Water Viscosity 
[1776-T 
u_(T) = |-—----— (Figure F-2) (2) 
Y. 2.65 (T) -89} 
Oil Viscosity 
" (fee ae 
ee) oie il E60) 6D (Figure F-3) (2) 


where the value of (v/) determines the sensitivity 


of the oil viscosity to temperature changes. 
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Kozeny-Carman Equation -- Its Application to 
Problems of Multiphase Flow in Porous Media", 
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(2)cGoTTFRIED, B.S.: "A Mathematical Model of Thermal Oil 
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APPENDIX G 


THE COMPUTER PROGRAM 
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THE COMPUTER .PROGRAM 


A computer program was written to solve the 


hot-water injection problem. The program may be divided 


into four eparts: 


(1) 


(2) 


(3) 


MAIN PROGRAM - The program will compute the 
pressure and saturation distribution. in a linear, 
One-dimensional, homogeneous reservoir model. 

This is accomplished by solving the continuity- 
Darcy flow equations for oil and water (Equations 
(11) and (12)) simultaneously, using the 
alternating direction explicit procedure (written 
for one-dimension) and the Newton-Raphson iterative 


technique. 


MATERIAL BALANCE SUBROUTINE - A cumulative and 
instantaneous material balance is carried out for 
each phase (oil and water) on a mass basis as 

a check on the main program results. Deviations 
from 1.0 will indicate the fractional error in 


the material balance. 


MATRIX FORMULATION FOR TEMPERATURE DISTRIBUTION - 
The pressure and saturation distributions obtained 


from the main program results are used to solve 
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the thermal energy equation (15). An equation 
must be written for every grid element since 

the equation is implicit. This subroutine forms 
the matrix of thermal energy equations (Crank- 
Nicolson Finite Difference Approximations) which 
must be solved to obtain the desired temperature 


distribution.at: the end of- the time interval. 


INVERSION OF A NON-SYMMETRIC MATRIX BY A MODIFIED 
SQUARE ROOT METHOD - This subroutine inverts the 
matrix which was developed by the subroutine 
discussed in Part (3). The temperature 
distribution is then calculated and used in the 
MAIN PROGRAM to calculate the temperature 
dependent properties for the next time interval. 
Mathematically, this subroutine performs the 
following calculations. If the thermal energy 
equations written for each grid element may be 


represented in matrix notation as, 


A T= b 


the. program will first. of all find the. inverse 


of A (assuming an inverse exists), then 
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COMPUTER PROGRAM .NOMENCLATURE 


CALL MATBAL 


CALL SOLUT = 


CALL TEMP = 


CPO (T) = 


CPOSTR=C* = 
Po 


CPW (T) = 


CPWSTR=C* = 
Py 


CPROCK=C = 
Des 
CUMMBO = 
CUMMBW = 
CW=C = 
W 
D8,D9,D12,D13 = 


Lo = 


call material balance subroutine 

call subroutine SOLUT to invert the 
matrix 

call the temperature subroutine to 
set up the matrix 

field unit conversion constant 

Oil compressibility 

function definition for specific heat 
of oil at constant pressure 


specific heat of oil at constant pressure 


function definition for specific heat 
of water at constant pressure 
specific heat of water at constant 
pressure 

specific heat of reservoir rock at 
constant pressure 

cumulative material balance of oil 
cumulative material balance of water 
water compressibility 

CONSTANTS 


same as LENGTH 


enistvoudue somsied Sabie gene: ; 
ous revel ot Toy onprngetpe Iiso 
| | xiadem 
oF snistuvordue e-tutesxoqmet aus ifss 
| xizzam odd qu tee. 


boy 


iy 


‘gnedenoo -notetevacs: Sav blot? 
ysibidtesesqmos Lio 

tse ottiosge x10t noitiniteb noitonut 
siwaestq tnstemoo ts Lio to 


eiueze1g tnstemoo ts Ifo to teen Sitioveqe . = 


teed oitiosgqe 10% aettitntves sotsonu® = 
e1weesta $nedenoo $s ustsw to. 
tastano> JI& esew 26 teen) ott ioege ia p It 
| orueestq 
$s Avot shay bets 20 tsad. shRhopae 
ewessxg: $nstenos io 
Lio to sonsisd isixedim ovitsiumss en 
t9t5w to eons ied tekved sm avisetume . 
ini xosew 


DKO (SO) 
DKW (SO) 
DT = At 


DX = Ax 


GP 


KE (SO) = 


KESTR 


KO(SO) = 


KW (SO) 


LENGTH 


MAXL 


MUO (T) 


MUOSTR 


MUW (T) 
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derivative of KO(SO) with respect to SO 
derivative of KW(SO) with respect to SO 
incremental time step size 
incremental space step size 
exponential 

mass injection rate of water 

number of grid elements 

function definition for effective 
thermal conductivity 

effective thermal conductivity at the 
initial oil saturation 

function definition for effective 
permeability to oil 

function definition for effective 
permeability to water 

cumulative space step size 

"dummy" length parameter to indicate 
reverse sweep (>L-Ax/2) 

function definition for oil viscosity 
viscosity of oil at initial reservoir 
temperature 


function definition for water viscosity 
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viscosity of water at initial reservoir 


temperature 


sensitivity of oil viscosity to temper- 


ature changes 

Oil accumulation 

mass Of oil injected = G_ At 

oil in place initially 
instantaneous oil material balance 
cumulative? oil* production (oil out) 
Oil produced during the nth time 
interval 

Oil produced during the (n+l) time 
interval 

absolute permeability 

hydrodynamic equation for the oil 
phase (function) 

derivative of PHI with respect to the 


iteration pressure PK 


derivative of PHI with respect to the 
iteration saturation SOK 


iteration pressure 
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PO = P. = initial pressure distribution 

POR = $6 = porosity 

PSTR = p* = reference pressure 

PZI = wv = hydrodynamic equation for the water 


phase (function) 

PZIP = Vp = derivative of PZI with respect to the 
iteration pressure PK 

PZISO = Vo = derivative of PZI with respect to the 


iteration saturation SOK 


ROO(P) = 0 = function definition for oil density 

ROOSTR = pe = density of oil at the reference 
pressure (PSTR = P*) 

ROROCK = aie reservoir rock density 

ROW(P) = ee function definition for water density 

ROWSTR = p* = density of water at the reference 
pressure (PSTR = P*) 

SOK = Sox = iteration saturation 

SOO = Soa = initial oil saturation 

SOXO = oil saturation at the inlet boundary (x=0) 

SWC = =e = connate water saturation 

SWO = Soi = initial water saturation 


TIME = real time 
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TLNGTH = L 
TO = Ty 
TXO = Ting 
WACCUM 

WIN 

WIPI 

WMB 

WOUT 


WOUTA 


WOUTB 


Y2 


YMAX 


= 150. - 


iteration temperature 

total length of model 

initial temperature distribution 
injection temperature at the inlet 
boundary (x=0) 

water accumulation 

mass of water injected = on At 

water in place initially 

instantaneous Water material balance 
cumulative water production (water out) 
water produced during the nth time 
interval 

water produced during the (n+l) time 
interval 

error in the guess of the iteration 
pressure PK 

error in the guess of the iteration 
saturation SOK 

maximum error allowed between two successive 
iterations in pressure before terminating 


iterations 
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maximum error allowed between two 
successive iterations in saturation 


before terminating iterations 
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SAMPLE COMPUTER INPUT - DATA CARDS 


(1) TO,S00,PO,GAS,SOX0,TXO 

(2) POR,ROROCK,CPROCK, PERM, SWO 

(3)  DT,DX,GP,D8,D9,CPOSTR 

(4) D12,D13,KESTR, ROOSTR,MUOSTR, ROWSTR 
(5) MAXL,TLNGTH,PSTR,SWC,CPWSTR,MUWSTR 
(6) NABLA 


(7) 8co7cw, CCl 


(8) NIN 

GL) . ¢p2.07-000 Dis) DO 14.700 Sais! 0.000 660.000 
(2) Of 40/1 65'..200 O2o0 \ LOOOn000 1.000 

(3) 0.010 Oro. 0 000.0) 1.000 0.000 0.650 
(4) iF 000 1.000 0.500 55.000 265.000 62.400 
C5) ih ele be 3.000 14.700 0.200 1.000 0.800 
(6) 0.560 


(7) S8800x10s"- 3520x1075: 2:637x107"% 


(8) 10 


DATA CHANGES: NABUAR=—0F 560,09) E20; sb 500 
DX = 0.150, 0.200, 0.300 


DT SOU PU COC) Ue Ua yr U6 ee 


et 
ig 


ATEWOA ATBOUM, reno aan 
aATaWUM, AT@WD , owe 8 


it oe f 
} Bids, a ms 
a , . ob : 
000.088 000.0 eis.0 ODT .84-~ O2TLO “000.088 N 


: 


000.£ 000,000L ES.0 o0s.205 08.0 
023.0 000.0  O00.f 000.02 defo” 090.0 
00.53 000,298 000.6% 002.0 000.£ \ ere 
008.0 000.4 008.0 Oor.eL © 000.€ ees 
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$ 
C 
C 
C 
C 
C 
C 
C 
C 
C 


IBFTC THACH 


LCC 
1Cl 
1C2 
104 
1c5 


1C6 
107 
108 
1Cc9 
110 
201 
2C2 


1 
* 
3 
4 
- 


#x*MATN 
AND SAT 
HOMOGEN 
SOLUTIO 
VARTABL 
CESCRIB 
ie ALS, 
CIMENSI 
EMPLCYE 
REAL MU 
FORMAT 
FORMAT 
FORMAT 
FORMAT ( 
FORMAT ( 
URE 
LK »65H 
HRS./ 
1H ,65H 
me cen mmf 
FORMAT ( 
FORMAT ({ 
FORMAT ( 
FORMAT ( 
FORMAT ( 
FORMAT ( 
FORMAT ( 
COMMCN 
COMMON 
COMMCN 
CCMMCN 
COMMON 
CUOMMCN 
***XxCATA 


REAC (5, 


REACT (5 
REAC (5 


REAC (5, 
REAC (5, 
Al eke 


REAC (5 
REAC (5 
#**xCCNS 
T2LCNG= 
XLNGTH= 
ILNGTH= 
IF (ABS 
JLNGTH= 
JT=JLNG 
COUT=C. 
wOUT=C 


wWINOUT 


DECK,NCLIST 
COMPUTER PROGRAM—-THE PROGRAM wILL CUMPUTE TRE PRESSURE 
URATION DISTRIBUTION IN A LINEAR, GNE-UIMENSIGNAL, 
ECUS RESERVOIR MODEL. THIS INVOLVES THE SIMULTANUDUS 
N OF TWO NON-LINEAR PARTIAL OLFFERENTIAL EQUATIUKS with 
—€ COEFFICIENTS. THE Tw CONTINUITY-DARCY FLOW EWUATIiiNS 
E THE SIMULTANEOUS FLGW OF TWO LIQUIDS, CIL AND WATER. 
ERNATING DIRECTION EXPLICIT PRCOCECURE (wiRkITTEN FOR UNE 
ON) AND THE NEWTGN-RAPHSON ITERATIVE TECHNIWUE wile 


D TO SOLVE THE PROBLEM. **#* 


O»MUWsKO eK Wy MUOSTR pNABLA,MUWSTReKROT »sLENGTHyMAXL yp KEdKESTX 


(6h. 4.3) 

(4E15.3) 

(318) 

Lh 928H MAXIMUM ITTERATICNS EXCECDED) 


1HK,y65H DISTANCE ~ PRESSURE WATER | — TEMPLRA 


TIME/ 
sa PSIA SATURATION K 


me ee eae ae 


) 

lh 9F8e4,F 14. 39F 15.69 Fl 6. 39F 1202) 

l1HJ,19H CIL PRODUCTION IS ,Fl2.4,8H POUNDS) 

TH »21H WATER PRGDUCTICN IS »,FlC.4,8H POUNLS) 

1H »40H CUMULATIVE MATERIAL BALANCE FGR OIL IS 4F13.4) 

lH »42H CUMULATIVE MATERIAL BALANCE FOR WATER Id oFLL.4) 
LK 943H INSTANTANEOUS MATERIAL BALANCE FUR JIL IS oF 1b.4) 
LH 945H INSTANTANEOUS MATERIAL BALANCE FUR WATER IS yt. 
TT(150)_SS(150) »PP( 150) pWENPL, TIME, JLNGT Hy DAgiT POR 
RCOSTR,»CU,sPSTRy»ROWSTReChyMUOSTR »pMUWSTR,SWCO,»,PERM,UGUT , WUT 
CUMMBO,CUMMbW,SCU,TLNGTHe PO,AIL(1L50) AS31150) pA 95y 95) 
AA(55,55) ,8(55) ,BB(55) AL (55,55) 7 AM(55,55) ,AU(55,55) 
AV(55955)9DIF(55)_,PCENT(55) 4 ¥(55) »NABLAY SW(150)9CCl 
OMB yWMBy COUTA,WOUTA,WINCUT,OCUIMB | 

READING *¥*#% 

100) TO,SCO,P0,GAS,SUXC,TXO 

2100) PORyROROCKyCPROCK yPERMy SWO 

2100) DT sDXyGPyD8,09,CPCSTR 

100) 012,01 3eKESTR»ROOSTR»MUOSTRyROWSTR 

1CO) MAXL»TLNGTHyPSTRySWCyCPWSTRyMUWSTK 

100) NABLA 

aL COUsCwW. CCT 

9102) NIN 

TANTS*#** 

TLNGTH-DX/2. 

TLNGTH/DX 

IF IX(XLNGTH) 

(XLNGTH-FLOAT(ILNGTH))-GE-0.9) ILNGTH=ILMNGTH] 
ILNGTH+I1 

Teel 

0 


0 
COUTA=0. 
WOUTA=0. 
=0.0 
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COUTMB=0.0 

CUMMBW=C0.C 

CUMMBEU=0.0 

CXDT=(POR¥*DX**2)/(CC1*DT) 

CCT=C.U 

XL=-C.5 

***PEYSTCAL PROPERTY EQUATIONS**¥* 
RCU(P)=RCOSTR¥EXP(CO*(P-PSTR) ) 
RCwW(P)=RCWSTR¥*E XP(CW*(P-PSTR) ) 
MUC(T)=MUCSTR*¥60. **NABLA* (T-46C. ) ®* (-NALELA) 
MUW(T)=(1776.-(T-460.))/(26.5*(T-460.)-89.) 
CPw(T)=CPhWSTR 
CPO(T)=CPOSTR*(C&+D9*%*(T-460U.) ) 
KE(SC)=KESTR*(DL24+D13*S0) 

Rod mG epr RM (1. —{ 1Lo-SO-SWC)/ 01 e-SWO)) 8¥2* (Le (01. -S0-SWO)/ 11 
lL SwC))**2) 
KwW(SC)=PERM*((1L.-SO-SWC)/(1.-ShC) ) *¥*4 

CKOA SG) =2 *PERMF(A Lem) --SG-SWC)/(L.-SwC) )*#2*(1.2-30-SHnCI/(1.- 
lp S Wee e2e(1.-—(01.-SO-SkC) /01Le-SWC)2*%2)*5S0/11.-5WC) #2) 
DKW(S0)=4.*PERM*((1.-SU-SWO)/(1L.-SWC) ) ¥*¥3*(-1L./41.-Sie)) 
***X FORWARD SWEEP*** 

*ee INITIAL CONDITIONS AND 6.Ce. AT K=0%*% 
NGP=IFIX(GP) 

CO 1 IT=lyNGP 

TT(1)=TO 

PP(I)=PO 

SS(1)=SOC 

SS(1)=SOXC 

TT(1)=TXCO 

SwW(1)=SWC 

wENP1=GAS 

*¥*¥*x BEGIN FURWARE SWeCP*** 

J=1 

ICUCE=0 

J=J+1 

COT=COT+1.6 

TIME=CT*¥LCCT 

SCNP1=SUXC 

TNPL=TXO 

PJN=PP(J) 

SCJN=SS(J) 

TIN=TT(J) 

PJPIN=PP(J+1) 

STIPIN=SSUJ+1) 

TIPIN=TT(U+1) 

KRCT=KO(SCNPL) 

IF (SCNP1.eLE-SOXC) KROT=0.6 

IF (ICODE-EQ.1) GOTOILILL 

XL=XL+1.C 

LENG TR=OX*XL 

IF (LENGTRK2GEeMAXL) XL=XL-1.0 

IF (LENGTFeGEeMAXL) J=J-1 

IF (LENGTR.GEeMAXL) GOTO9 

GCVC1L1L2 

XL=XL-1.C 
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LENGTR=DX*XL 
L112) AJN=KO(SCJIN) /MUG( TIN) 

HBIPIN=KO(SOJPIN) /MUO(TJPIN) 
PIN=KW(SCIN) /MUW( TIN) 
BUPIN=KW(SUJPIN) /MUW(TJPLN) 
ANPL=KG(SONPL)/MUO(TNPL) 
BNPL=KW(SONPL)/MUW(CTNPL) 

C ***TTERATICN FOR PK AND SUK*** 
K=0 
TK=TJN 
SOK=SS(J)-U.9TL¥*DT 
PK=PJIN¥*1.2 

3 K=K+1 

IF (K.GT.20) GUTO5 
UL=RCOUPK ) 
U2=KO(SCK) 
IF (SCK.LE.SUX0) U2=C.0 
IF (SCK.CE.0.8) U2=PERKM 
U3=MLC(TK) 
L4=RCW(PK ) 
U5=KwW(SCK) 
IF (SOK.CE.0.8) U5=C.0 
[F (SOK.LE.SUXU) US=PERM 
U6=MUW(TK) 
DU2SCK=DKO(SCK) 
IF (SCKe.LE~SOXO0) DU2S0K=0.0 
IF (SCK.GE.U0.8) OU2SOK=uU.0 
CUSSCK=DKW( SOCK) 
IF (SCK.GE.0.8) OUSSOHK=06.0 
IF (SCK.LE~SCX0O) SUK=SOX9 
IF (J-GE.2) R1=1.0 
IF (J.G5E22) R2=1.0 
IF (LENGTR.GE-T2LONG) RL=4./3. 
IF (LENGTEF.CE.T2LONG) R2=8./3. 
IF (J.EQ.2.-AND.ICODE.EuL1) GUTCLY 
IF (J.E@2.2-AND~-ICCDE.EQ.0O) GOTLI8 
IF (ICODE.«EUeL) GOTO16 
IF (JeGT.22) CELTAP=(PNPI1-PK) 

C **eENEWTON-KAPHSON TTERATIVE TECHNI wut FCR INTERTUR GRIGG PUINIS 
(J GREATER THAN 2) AND FOR THE OUTLET BOUNDARY UN Thi FCRWARKI, 
SWEEPS SOM * * 

PHI=R1*¥(L2/U34ANP1) *DELTAP/2.+R2*(AINtASPIN) ¥(PSPLWRPINI J 2.-1 Ki * 

L (SGK-SOJN+CO*( SOK+SOJN) *(PK-PJN)/2.). 

PZI=ERL*¥(U5/U64BNP1) *¥DELTAP/2.4+R2¥*(BINtK SPIN) ¥(PIPLN-PU WIS 2 e-t XT 

lL (SCJN-SCK+CW*(1L.-CSOK+SOJSN)/2.2)¥*(PK-PIUN)) 

C *eXELEMENTS OF ITERATICN MATRIX*** 
PRIP=RL¥*(L2/U3+ANPL)/ (0-2. )-DALT*(CO*¥(SUK+SUUNI/ 2.) 
PZIP=RL*¥(LS/U64BNPL)/(-2. )-DADT*¥(Cw*#(L.-(SUK+SOUNI/ 2.) ) 

24 PRISC=R1*0U2SCK*DELTAP/(2.*U3)-CXDT*(1L.+C0%*(PK-PINIJ ec.) 
PZISC=R1L*DUSSOK*DELTAP/(2.¥*U0) tOXDT*¥( Let Cw*( PK-PIN)I/ 2.) 
GOTCSC 
*&*NEWTUN-RAPHSON ITERATIVE TECHNIQUE FOR THE FNLET GUUNDARY 
(J EQUAL TO 2) ON THE FORWARU SWEEP. *¥** 

18 PRHI=(AJIN+AJPIN) *(PIPLN-PIN) /2.-DXDI*( SUK-SUINFCU* (SOK + SUNN ) * 

L (PK-PJUN)/2.) 
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PZI=wFNP1L¥*¥DX/(CCL¥ROW(PK) )+(BIJNtBUPIN) *(PJPLN-PJN)/2.-DXDT* 
(SOJN-SCK4+CW*(1.-(SOK+SOJN)/2.) *(PK-—PJN) ) 

PHIP=-DXCT*(CO*(SOK+SOJN) /2.) 

PZIP=-DXDT*(CW*¥(1.-(SOK+SOJN) /2.)) 

PHISC=-DXOT*(1.+CO*(PK-PJN)/2.) 

PZISO=DXCT*(1.+CW*( PK-PJN)/2.) 

GCTCSO 

***NEWTON-RAPHSON ITERATIVE TECHNIQUE FOR THE INLET &BUUNUARY 

(J EQUAL TO 2) ON THE REVERSE SWEEP. *%* 

PHT=(U2/U34+AJPIN)*(PJIPLN—PK)/2.-DXDT* (SOK-SOJUN#CU* (SUK +SGJIND % 
(PK-PJN)/2.) 

PZL=WENP1L*DX/(CC1L¥ROW( PIN) )+(U5/U64B5NPLN) *(PIPLN-PK)/2.-DXDT * 
(SCJN-SCK4+CW*( 1L.-(SOK+SOJUN)/2.)*(PK-PJN) ) 

GOTCGS98 

***NEWTON-RAPHSON ITERATIVE TECHNIQUE FOR INTERIGR GRID POINTS 
(J GREATER THAN 2) AND FOR THE OUTLET BOUNDARY GN The REVERSe 

SWEEP. ¥#%% 

IF (J.GT.22) DEL TAP=(PNPL-PJN) 

PHI=R1L*(AJN+ANP1)¥*¥OELTAP/2.+R2*(U2/U34AJPLN) ¥(PIPIN-PK)/20-UXUI * 
(SCK-SOJN+CO*( SOK+SOJN) *(PK-PUN) /2.) 

PZIT=R1L*(BIN+BNPL) *DELTAP/2.4R2¥*(US/U64BUPIN) *(PUPIN-PK)/20.-LKUT* 


L (SCJN-SCK+CW*( 1.—(SOK+SOUN)/2-) *(PK-PJN) ) 


#*KELEMENTS OF ILTERATION MATRIX¥#* 

PHIP=R2*(U2/U3+AJPLN) /(-2.)-DXCT*(CO* (SGK+SUSNI/26) 
PZIP=R2*(US/U6+K SPIN) /(-2.)-OX0T*(CW*(1L.-(SUK#SUINI/2.)) 
PHI SO=R2*0U2 SOK * (PUP IN-PK)/(2.*U3)-UXDI*( 1. +CO*(PK-PJNI/ 26) 
PZ1SC=R2*DU5SOK* (PUP IN-PK)/ (2. #U6) DXDT *( 1. +CW#(PK-PINI/ 26) 
YL=(PZI*PHISO-PHI*PZISO)/(PZISC#PHI P-PZIP*PHISU) 
¥2=(PZL*PHKIP-PHI*PZIP)/(PZIP*PHISO-PHIP*PZISU) 

YMAX=1.0E-04 

YMAX1=1.CE-06 

IF (ABS(Y1/PK).LEsYMAXsANDeABS(Y2/SOK).LEsYMAXL) GOTL2 
PK=PK4Y1 

SOK=SOK+Y2 

GOTC3 

CONTINUE 

IF (J.GT.2) GOTC187 
PP(1)=DELTAP+PK 
PR=PP(1) 

PP(J)=PK 

SS(J)=SOK 
SWK=1.-SCK 
SW(J)=SWK 

IT(J)=TK 
RCCP=(U4¥*CPW(TK) *(L.—-SOK) +UL*CPU( TK) SUK) *#PUR+ROROCK*CPRUCK*(1L.- 


Ll PGR) 


ACV=(UL*¥CPO(TK) *U2/U3+U4*CPW(TK) ¥U5/U6) 
ASL (J) =26#ROCP*DX**2/(CCL*ACV*CT) 
AJS3(J)=KE(SGK)/(CCL*ACV) 

IF (ICODE.EQ.1) J=J-1 

JS=J 

IF (JeEQ.1.AND.ICOQDE.EU.1) XL=XL-1.0 

IF (Js£Q.1.AND.ICODE~EW.1) GOTC185 

IF (ICODE.EQ.1) GOTI93 

PNP 1=PP(J) 
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SONP1=SS(J) 

TNPL=TT(J) 

J=J4l1 

GOTQ4 pe 
CALL MATPAL 

CALL TEMP 

CALL. SOLUT 

NCDT=IFIX(CDT) 
NCDT2=NDCT/NIN 

IF (NCNT2.EQ.1) GOTO18O 
J=JS 

GOTO20 

WRITE (6,105) 

DIST=0.0COO 

CO. 4C0. Jw=l,JT 

WRITE (6,106) DIST,PP( JW) »SWwldw) pT (JW) eo TIME 
XX=EX 

TPitelWe le Lele IWec Wel JI—1) 1, XXEDa/ Ce 
DOIST=CIST+Xx 
CONTINUE 

WRITE (6,107) OOQUT 
WRITE (6,108) WOUT 
WRITE (6%109) CUMMBO 
WRITE (6,110) CUMMBW 
WRITE (6,201) OMB 
WRITE (6,202) WMB 
NIN=NIN#I1C 

J=J5 

GCTO2C 

€*XREVERSE SWEEP *** 
ICOrEs, 

COT=CCT+1.0 
TIME=CT*COT 
{teed eS OTe 

KIK=J+4l 
SS(KIK)=SS(¢J) 
Sw(KIK)=Le-SS(KIK) 
PNPL=PP(J-1) 
SCNPL=SS(J-1) 
TNPL=TT(J-1) 

GOTC4 

WRITE (6,104) 

CALL EXIT 

ENC 
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SUBRCUTINE MATBAL 

***MATERTAL BALANCE SUBROUTINE--A CUMULATIVE AND INSTANTANEOUS 
MATERIAL BALANCE IS CARRIED OUT FOR EACH PHASE (UTIL ANU wATE?) 
AS A CHECK CN THE MAIN PRCGRAM RESULTS. DEVIATIONS FRUM 1.0 wiILl 
INCICATE THE FRACTIONAL ERROR IN THE MATERIAL BALANCE. *** 
REAL MUQsMUWsKO eKWe MUOSTR yNABLAy MUWSTRyKRGT »yLENGTHpMAXL eg KE gK ESI 
COMMON TT(150)_,SS(150) »PP( 150) sWENPL> TIME, JLNGTH,DX90T y POR 
COMMCN RCCSTRsCOsPSTR»RGWSTReCWweMUUSTRyMUWSTRe Sly PERM GOUT pwlul 
COMMCN CUMMBO »CUMMBW SOG» TLNGTHe POs AJL (150) pAU3( 150) 9A(559 95) 
COMMON AA(55,55) 4B(55) »BB (55) 9 AL (55955) 9 AM(55955) ,AU(55y 95) 
COMMCN AV(55,55) eDIF(55),PCENT(55),Y(55) sNABLA,SW(150) ,COl 
COMMCN OMByWMByOOUTA,WUUTA,WINCUT »OGUTMB 

FORMAT (lt ,8E15.5) 

RCO(P)=RCCSTR*EXP(CC*(P-PSTR) ) 

ROW(P)=ROWSTR¥EXP(CW*(P-PSTR) ) 

MUG (T) =MUCSTR*60.**NABGLA* (T-46C.~ ) ** (-NABLA) 
MUW(T)=(1776.-(T-460.) )/(26.5%*(T-460. )-89.) 

KC (SO) =PERM*(1.-(1.-SC-SWC)/(1.-SWC) ) ¥*®2*(1.-(01.-SU-SwWC)/(1.- 
1 SwC))**2) 

KWw(SC)=PERM*((L.-SO-SWC)/(1--ShC) ) **4 

OIN=0.0 : 

CIPI=TLNGTH*¥POR*ROO (PO) *SOG 

wIN=hFENPL*TIME 

WIPIL=TLNGTH*®POR*ROW (PO) *(1.-SCC) 

wACCUM=0.C 

CACCUM=0.C 

CO 1SC I=2,JLNGTH 

Pes hl) 

SiaSS(E) 

PL=PP(1) 

QACCLEM=OACCUM4D X¥POR¥*SL¥RCO (PL) 

wWACCUM=WACCUM4D X¥POR¥*(1.2-SL)*RCWI(PL) 

PLM1=PP(JLNGTH) 

PL=PC 

PR=(PLML+PL)/2. 
wWOUTP=(RCW(PR)*KW(ISL) /MUW(TL) ) *¥(PLMIL—-PL) ¥2.*DT*CC1/(UX/2.) 
COUTB=(RCO(PRI*KO(SL)/MUU(TL) ) *(PLML-PL)*2.*DT*CCI1/(0X/2.) 
WOLTAV=(hCUTA+WCUTB) /2. 

COLTAV=(CCUTA+OCUTB) /2. 

IF (TIME.EQ.2.*0T) WOUTAV=WOUTE 

IF (TIME.EQ.2.*NT) GOUTAV=O0UTB 

IF (TIME«EQe2e¢*0DT) HACCUM=WIFI 

IF (TIMEEQ.2-*CT) GACCUM=CIPI 
wMB=(WACCUM-HACCUM) / (WENP1¥*2.*CT-WOUTAV) 
CMB=(GACCUM-OACCUM) /COUTAV 

FKACCLUM=WACCUM 

GACCLM=OACCUM 

wOUT=wOUT+WOLTAV 

COUT=COUT+COLTAV 

COUTMB=OIPI-GACCUM 

wINCLUT=WwACCUM-WIPI 

CUMMBW=WINOQUT/( WIN-WOUT ) 

WOUTA=WOLTE 

COUTA=O00LTB 
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IF (QOUT.LE.0.0) GOTC186 
CUMMBC=OCLTMB/OOUT 
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¢IBFTC MAYRE CECKysNCLIST 
SUBRCLTINE TEMP 
#e*MATRIX FORMULATICN FOR TEMPERATURE OISTRIBUTIGN-- THIS 
SUBRCUTINE USES THE SATURATION AND PRESSURE DISTRIBUTIGNS 
OBTAINED FROM THE MAIN PROGRAM TG FURM A MATRIX UF THERMAL 
ENERGY EQUATIONS (CRANK-=NICULSCN FINITC DIFFERENCE EWUATIONS) 
WHICt WHEN SOLVED WILL YIELD A NEW TEMPERATURE DISTRIBUTIGN. *** 
REAL MUO »MUWsKO KW, MUOSTRyNABLAyMUWSTReKROI pLENGTH ey MAAL KE ph ST? 
CUOMMCN TT(150),S5S(150) »PP(150) pWENPL» TIME sJLNGTHyOAgiil yPuR 
COMPCN RCOSTReCCyPSTReROWSTRyChyMUUSTReMUWSTRy SW PERM SCOUT yp viJul 
COMMEN CUMMEO yCUMMBW SCO, TLNGTHyPO,AJL(1L50) pAS3(15U) 74199955) 
COMMON AA(55455) 48155) 986 (55) 9 4L (55455) ¥AM(55455) yAU(5 9555) 
COMPMCN AV(55455),0IF (55) »PCENT(55) 9 ¥(55) yNABLA,SHIL5U) »CUl 
COMMCN OMB,WMByCOUTA,WOUTAyWINCUT »GQUOUTMB 
1C1 FORMAT (1h »8E15.5) 
JJLNG=JLNGTH+1 
CO 3S I=l»yJJLNG 
CG 39 J=lyJJLNG 
3S A(1I,¥)=O08C0 
MLNGTH=JLNGTH-1 
322 
I=1 
7S AJS2=-PP(J42)+6.¥*PP(J4+1)-—3.*PP(J)-2.¥*PP(J-1) 
AJS4=-TT( 542) 46.¥*TT(J41)-36*TT(J)-2.*TT(J-1) 
AJSS=TT(S41)—-2.*TT(J) 4TT(U-1L) 
BBL=AJS1(5)4AS2/12.4+2.¥*AJS3 (J) 
BB2=-AJ2/6.-AJ3 (J) 
BB3=AJ2/36. 
BB4=AJ2/18.-AJ3 (J) 
RBBS=AJL(JS)*TT(J)+BB3*AJS44AS3(5)*AJ5 
IF (JeEQ.2) GOTO78 
IF (JeEQ.3) GOTOTT 
IF (J.EQ.4) GOTUT6 
GCTG7S5 
LEL AS2=20¥*PP(J41)4+3.¥*PP(J)-66*PP(J-1)+PP(J-2) 
ASG=2.¥*TT(J41)¢30%TT (J) -60¥*TT(J-L) + FT (J-2) 
AJS5S=TT(JS41)-2.* TTI) +TT(U-1) 
PB2=AJ1L(JS)—-AJ2/12.4+2.¥*AI3 (J) 
BB3=-AJ2/18.-AJ3( J) 
BB4=-AJ2/36. 
PRL=AJ2/64-AJ3(J5) 
BBS=AJL( J) *T T(J )-BB3*AJ44A53 (5) *AN5 
GCTC140 
7@ A(l,1)=BBl 
A(l,2)=B8B2 
A(1l,3)=BB3 
B( 1) =BB5-BB4*TT(J-1) 
J=J+1 
l=1+1 
GCTC79S 
77 A(2,1)=BB4 
A(2,2)=BB1 
Al2,2)=BB2 
A(294)=RB3 
B(1)=8B5 
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J=J+1 

l=I +1 

GOTC7S 

L=3 

MMM=C 

M=] 

NN=1 

M=M4NN 

MM=M—-MMM 

IF (MM.LEQ.2) GOTN43 
IF (MM.EC.3) GOTO45 
IF (MM.EC.4) GOTO46 
IF (MM.ZEQ.5) GOTO47 
GCT048 


A(lL»M™)=BB4 
GOTC44 
A(L,»M)=BBl 
GCTC44 
A(t »M)=BR2 


IF (JeGEeMLNGTH) Al(L»M)=8B2+6B3 
IF (JeEQeMLNGTH) Bt I)=8B5 


IF (JeEQ.MLNGTH) GOTO42 
IF (Je.EQeJLNGTH) GOTO74% 
GOTC44 

A(L»M)=BB3 

B(I)=bB5 

GCOTC44 

J=J+1 

I=I +1 

IF (J.EQeULNGTH) GOTOL61L 
GoTC79 

L=L+1 

NN=L-2 

M=1] 

MMM=MMM41 

GCOTC49 

L=L+1 

NN=L-3 

M=1 

GOTC4S 

RETURN 
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SUBROUTINE SOLUT 

*** INVERSION OF A NON-SYMMETRIC MATRIX BY A MODIFIED SQUARE FUuI 
METHOC--THIS SUBROUTINE INVERTS THE MATRIX SET UP BY SUBROUTING 
TEMP AND THEN SOLVES FOR THE NEW TEMPERATURE DISTRIBUTION 

wHICH IS THEN USED IN THE MAIN PROGRAM TO CALCULATE THE TEMPERATURE 
CEPENCENT PROPERTIES FOR THE NEXT TIME INTERVAL. ¥*¥** 

REAL MUO »MUWsKOeKWe MUOSTR pNABLAy MUWSTRyKROT pLENGTHyMAXL ep KEsKi S12 
COMMON TT(150)9SS(150) »PP(15C) sWENPle TIME sy JLNGTHyDAyCT pPOR 


COMMON RCOSTRyCOyPSTRyROWSTReChyMUUSTRyMUWSTRe SWC yPERM OOUT pwCUl 
COMMCN CUMMBCeCUMMBW, SCUyTLNGTHyPO,sAJLI1L50) AS30150) 9A( 55955) 
CCMMON AA(55,55)9B(55) BB (55) 4AL (55455) 9AM(55955) yAU( 59455) 
COMMON AV(55,55) »DIF (55) yPCENT(55),Y(55) pNABLA,SW(15C) CCL 
COMPMCN OMB,WMByCOUTA,WOUTA,WINCUT »OOUTMB 
FORMAT (1F ,»8E15.5) 

FORMAT (14X,14HSOLUTION ERROR/14X,6F15.6) 
CALL FPTSW 

M=JLNGTH-L 

CC l I=1,W 

CO 2 J=lyé 

AAC Ts JI=HA(1, JS) 

eEB(1T)=6(1) 

CALCULATICN OF L AND U 
COelCSxK=2,é 

KK=K-1 

CDC 1C4 J=1,KK 

ALtdsk)=C.0 

AM(J9K)=020 

AU(K,J)=C.0 

AV(K,J)=C.0 

CONTINUE 

IF(A(1l,1).EG@-0-) GO TO 150 
XIS=ABS(A(1,1L)) 
AL(1l,1)=SQRT(XIS) 
AU(LeLI=HAClLelI/SAL( 191) 
CG 10? J=2_9™ 
AU(L,»J)=A(L,IISAL( 151) 
AL( Js LIHAAC Je, LI/SAU(1L»1) 
COO ALS  K=2eM 

KK=K-1 

VALUE=0.C 

CN 108 J=1l,KK 
VALUE=VALLE+AL(K_J) *AU(J9K) 
Z2X=A(K»yK)-VALUE 
IF(ZXe-EQ.Ce) GO TO 150 
ZXX=ABS(ZX) 
AL(K»eK)=SQRT(ZXX) 
AU(KyK)=ZX/AL(K 9K) 

KP=K+1 

KK=K-1 

IF(KP.GT.M ) GO TO 115 

CO 112 1=KP eM 

Z¥2#Ce0 

Zw=C.0 

CO 110 LP=1l,KK 
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ZV=ZV+AL(KyLP)*AU(LP 91) 
LLC ZW=ZWHAL( 1 yLP)*AU(LPoK) 
AU(KyT)=(A(Ky1)-ZV) /AL (Ky K) 


112 AL(IsK)=(AC I yK)-—ZW) /AU(K,K) 
115 CONTINUE 
L ANC U ARE CALCULATED 
PROCEEDING TO CALCULATE L-INVERSE AND U-INVERSE 
DO 119 K=lyM 
AM(KyK)=1.0/AL{( KK) 
11S AV(KeK)=1.0/AU( Ko K) 
CO 125 K=2 9M 
KK=K-1 
Giel220215KK 
ZQ=C.0 
DG 120 L=JsKK_ 
12C ZQ=ZG+AL(KyL)#AM(LyJ) 
122 AM(KyJ)=-ZQ/AL( KyK) 
12& CONTINUE 
IJN=M+1 
CO 135 KL=lyIJN 
BS UNAKL oo Mee 
IF(KeLE.1) GO TCO 136 
126 DO 130 JK=l»K 
J=K-JK 
PFAU.T.1) 60-101135 
127 ZR=0.0 
JP=J+1 
DO 128 L=JP,yK 
12& ZR=ZR+AU( Jol) #AVILYK) 
13C AV(JsK)=-ZR/AUL JJ) 
135 CCNTINUE 
PROCEEDING TO CALCULATE A-INVERSE=U-INVERSE*L-INVERSE 
136 DO 140 K=l eM 
DC 140 J=lyM 
A(JK)=0.0 
CO 137 L=1l,™ 
A(SeK)=Al(J9K) tAV(I Sel) ¥AM(L 9K) 
137 CONTINUE 
14C CONTINUE. 
GO Tc 152 
15C WRITE(6,151) 
151 FORMAT(1E »29H THE METHOD IS NCT APPLICABLE) 
152 CONTINUE 
DO 153 I=lyM 
¥(1)=C. 
CO 153 J=leM 
153. Y(IT)=Y(1)4AC1T,J)*B(J) 
CO 154 I=1,JLNGTH 
KK=[+1 
154 TT(KK)=Y(1) 
JIJ=JLNGTH+1 
TT(JIJ)=TT( JLNGTH) 
CALL FPTSW 
CO 9 J=lyM 
BOB=C. 
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CO 12 I=1l,™ 
BOB=BCB+AA(Je1)¥*Y(1) 


CIF(J)=BB(J)-808 


PCENT(J)=DIF (J) *100./BB(J 
CONTINUE 

RETURN 

END 
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COMPUTER “OUTPUT I 
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NABLA 0.560 


Ax = DX 
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DISTANCE 
FEET 


_—— ee we ws oe owe 


2-8500 
3.0000 


PRESSURE 
PSTA 


ee ee ee ee eee a 


25-448 
25.066 
24-686 
23.791 
22.537 
21.330 
20-124 
18.919 
17.713 
16.508 
15.303 
14.700 


OIL PRODUCTION IS 
WATER PRODUCTION IS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 


INSTANTANEOUS MATERIAL BALANCE FOR OIL 
INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 


OISTANCE | 
FEET 


0.COCC 


PRESSURE 


PSTA 


2-850C 
3-COCC 


GIL PRCCUCTICN IS 
WATER PROCUCTION IS 


CUMLLATIVE MATERIAL BALANCE FOR OIL IS © 


CUMLLATIVE MATERIAL BALANCE FOR WATER IS 


INSTANTANECUS MATERIAL BALANCE FOR OIL IS 
INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 


WATER TEMPERATURE 
SATURATION R 
1.000000 520.000 
0.482922 520.000 
0.407252 520-000 
0.294014 520.000 
0.251132 520.000 
0.250006 520.000 
0.249998 520.000 
0.249999 520.000 
0.250000 520.000 
0.250001 520.000 
0.250001 520.000 
0.250001 520.000 
1.6944 POUNDS 
0.0078 POUNDS 
1.0146 
1.0010 
IS 0.9995 
1.0051 
WATER  ————_-§- TEMPERATURE 
SATURATION R 
1-CCOCOO 660.000 
0.464256 589.494 
__0-407581 aici Dee 
0. 303895 527-804 
0.251841 521.927 
0.25C015 520.391 
0.245999 520.061 
0.249999 520.001 
__0.25C000 519.992 
~ 0.259001 519.994 
0.250002 519.999 
0.25C002 519.999 
1.6328 POUNDS 
0.C075 POUNDS 
E FOR ¢ ~ 1.0343 
0.9830 
1.0112 
0.9682 


TIME 
HRS. 
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DISTANCE PRESSURE WATER TEMPERATURE 


FEET PSIA SATURATION R 
0.0000 21.992 1.000000 520.000 
0.1500 21.815 0.525016 520.000 
0.4500 21.608 0.474276 520.000 
0.7500 21.299 0.437938 520.000 
1.0500 20.795 0.355678 520.000 
1.3500 19.881 0. 262356 520.000 
1.6500 18.928 0.250148 520.000 
1.9500 17.989 0.250003 520.000 
2-2500 17.049 0.250002 520.000 
2.5500 16.110 0.250003 520.000 
2-8500 15.170 0.250004 520.000 
32-0000 14.700 0.250004 520.000 

OIL PRODUCTION IS 321353 POUNDS 
WATER PRODUCTION IS 0.0144 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0163 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.9256 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 1.0222 


INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.8016 


DISTANCE PRESSURE WATER TEMPERATURE 
FEET PSIA SATURATION R 
0.COoCcOo | 22.132 1.660060 ~~ 660.000 
0.1500 22.042 0.5C0967 609.092 
0.4500 21.93C 0.463910 567.487 
0.7500 21.741 0.442757 541.779 
1.C500 21.361 0.389864 528.778 
(1.3500 20.535 0.284928 523.159 
1.6500 19.458 0.250794 520.999 
1-$500 18.404 0.250011 520.271 
222500 17.347 0.250003 520.058 
2.5500 16.288 0.250003 520.009 
2.8500 15.23C 0.250004 520.012 
3.C000 14.70C 0.250004 520.012 
CIL PRCCUCTICN IS 322311 POUNDS 
WATER PROCUCTICN IS 0.0148 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0201 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.9577 
INSTANTANECUS MATERIAL BALANCE FOR OIL IS 1.0048 


INSTANTANECUS MATERIAL BALANCE FOR WATER IS 0.8880. 
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DISTANCE ———— PRESSURE WATER TEMPERATURE 
FEET ji me 1-4 Cons SATURATION ~ R 
0.0000 20.791 1.000000 520.000 
0.1500 20.663 0.548802 520.000 
0.4500 20.504 0.494805 520.000 
0.7500 20.291 0.469797 520.000 
1.0500 20.019 0.450668 520.000 
1.3500 19.606 0.397898 520.000 
1.6500 18.869 0.291054 520.000 
1.9500 17.915 0.251095 520.000 
2.2500 16.995 0.250014 520.000 
2.5500 16.077 0.250005 520.000 
2.8500 15.159 0.250006 520.000 
3.0000 14.700 0.250006 520.000 
OIL PRODUCTION IS 4.4804 POUNDS 
WATER PRODUCTION IS 0.0205 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0195 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.8841 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 1.0258 


INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.7855 


DISTANCE PRESSURE WATER TEMPERATURE 
0.COCO 22-5CC 1.0C0C00 660.000 
0.15C0 22.384 0.521968 618.755 
0.45CO 22.306 0.484723 581.751 
0.7500 22. O8G 0.467383 554.304 
1.0500 21.775 0.455531 537.079 
1.3500 21.368 0.435051 527.752 
1.6500 20.703 0.364465 523.242 
1.550C 15.458 0.266967 521.249 
2.2500 18.084 0.250245 520.434 
2.5500 16.730 0.250008 520.139 
2.85C0 15.374 0.250006 520.038 
3.C000 14.700 0.250006 520.038 

CIL PRCOUCTICN IS 4.8443 POUNDS 

WATER PROCUCTION IS 0.0222 POUNDS 

CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0181 

CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.9550 

INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 1.0180 


INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 1.1458 
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DISTANCE PRESSURE WATER 
FEET PSIA SATURATION 
0.0000 19.594 1.C00000 
0.1500 19.490 0.566115 
0.4500 19.357 0.509379 
0.7500 19.182 0.481683 
1.0500 18.965 0.467172 
1.3500 18.714 0.456064 
1.6500 18.381 0.423168 
1.9500 17.802 0.328627 
222500 16.899 0.255227 
2.5500 16.016 0.250059 
2.8500 15.139 0.250008 
32-0000 14.700 0.250008 


OIL PRODUCTION IS 527694 POUNDS 
WATER PRODUCTION IS 0.0265 POUNDS 
“CUMULATIVE “MATERIAL BALANCE FOR OIL TS 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 
INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 
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1.0206 
0.28545 
1.0236 
0.7488 
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_CISTANCE PRESSURE WATER * 
FEET PSIA SATURATION 
0.COCC 15.594 1.-C00C00 
0.1500 15.955 0.536567 
_ 0.4500 15.893 0.501980 
0.75CC 15.821 0.481531 
1.C500 15.75S 0.469052 
1.3500 15.694 0.460466 
1.6500 15.619 0.451330 
1.95CO 15.514 0.417384 
2-25C0 15.332 0.318794 > 
2-5500 15.077 0.253643 
2.8500 14.831 0.250058 
3.COCO 14.70C 0.250058 
CIL PRCCUCTICN IS 6.3702 POUNDS 
WATER PRODUCTION IS 0.C292 POUNDS 


CUMULATIVE MATERIAL BALANCE FOR WATER IS 


INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 
INSTANTANECUS MATERIAL BALANCE FOR WATER IS 
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DISTANCE PRESSURE WATER TEMPERATURE TIME 
FEET PSIA SATURATION R HRS. 
0.0000 18.583 1.C00000 520.000 45.00 
0.1500 18.494 0.579732 520.000 45.00 
0.4500 18.379 0.521341 520.000 45.00 
0.7500 18.227 0.491833 520-000 45.00 
1.0500 18.037 0.474753 520.000 45.00 
1.3500 17.817 0.465370 520.000 45.00 
1.6500 17.575 0.458321 520.000 45.00 
1.9500 17.283 0.437631 520.000 45.00 
2.2500 16.810 0.363214 520.000 45.00 
2.5500 15.999 0.265950 520.000 45.00 
2.8500 15.128 0.250295 520.000 45.00 
3.0000 14.700 0.250292 520.000 45.00 
OIL PRODUCTION IS 720126 POUNDS 
WATER PRODUCTION IS 0.0322 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0210 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.8310 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 1.0226 
INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.7284 
CISTANCE PRESSURE? i WATERS. 5 TEMPERATURE Liaise 
FEET PSIA SATURATION R HRS. 
0.COCC 13.218 1.-CCCCOO 660.000 45-200 
0.15CC 13.216 0.548364 628.759 45.00 
O24500 i “(ssi TI 0.504980 Sy sche. t,o aoe 
0.7500 13.16C 0.481586 571.448 45.00 
1.C5CO KIST97 0.468247 551.617 45.200 
1.35CC 13.265 0.459386 538.095 45.00 
1.65CCO 13.359 0.448716 529.678 45.00 
1.9500 12.454 0.4C6644 524.856 45.00 
222500 V2s74S (“AIT TST CC aaa ee 45.200 
2-55CC 14.135 0.251512 yi, | y ae Pees 
2-85CC 14.513 0.250035 518.668 45.00 
3.C000 14.700 0.250035 518.668 45.00 
CIL PRCCUCTICN IS 622C22 POUNDS 
WATER PROCUCTION IS 0.C284 PGUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS Naren Ne ey 
CUMLLATIVE MATERIAL BALANCE FOR WATER IS 0.7428 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 0.8477 
INSTANTANECUS MATERIAL BALANCE FOR WATER IS -0.2606 
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COMPUTERY OUTPUT, I1 


Y = NABLA = 1.120 
praca? ABD. = 0.300 feet 
Meo Ee ay = 0.020 Hrs. 
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DISTANCE PRESSURE WATER TEMPERATURE TIME 
FEET PSIA SATURATION R HRS. 
0.0000 25.341 1.000000 520.000 9.0 
0.1500 24.965 0.482909 520.000 9.0 
0.4500 24.585 0.406842 520-000 9.0 
0.7500 23.706 0.293873 520.000 9.0 
1.0500 22-463 0.251126 520.000 9.0 
1.3500 rit Wg oo | 0.250005 520.000 9.0 
1.6500 20.073 0.249997 520.000 9.0 
1.9500 18.879 0.249998 520.000 9.0 
222500 17.685 0.249999 520.000 9.0 
2.5500 16.491 0.250000 520.000 9.0 
2.8500 15.297 0.250000 520.000 9.0 
3.0000 14.700 0.250000 520.000 9.0 
OfL PRODUCTION IS 1.6768 POUNDS 
WATER PRODUCTION IS 0.0077 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0239 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.9996 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 0.9962 


INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.9930 
DISTANCE PRESSURE WATER TEMPERATURE TIME 
FEET PSIA SATURATION R HRS. 
0.0000 23.822 1.000000 660.000 9.0 
0.1500 23-636 0.474959 589.324 9.0 
0.4500 232437 0.401498 545.928 9.0 
0.7500 22.805 0.291652 527.817 9.0 
1.0500 21.753 0.250992 521-930 9.0 
3500 20.682 0.250004 520.394 9.0 
00 19.597 0.249998 520.065 9.0 
1500 18.509 0.249999 520.005 9.0 
-500 17.421 0.249999 519.996 9.0 
225500 16.332 0.250000 519.996 9.0 
28500 15.2244 0.25C000 519.999 9.0 
3.0000 14.700 0.250000 519.999 9.0 
GIt PROOUCTION IS 1.5253 POUNDS 
WATER PRODUCTION IS 0.0070 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0853 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.9633 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 1.0307 
INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.9353 
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DISTANCE PRESSURE WATER TEMPERATURE TIME 


FEET PSIA SATURATION R HRS. 
0.0000 22.889 1.C00000 520.000 18.0 
0.1500 22.718 0.524720 520.000 18.0 
(0.4500 —_ 22.511 0.474521 520.000 18.0 
057500 ~~~ ~~ 22.226 0.444946 520.000 18.0 
1.0500 21.730 0.379558 520.000 18.0 
1.3500 20.764 0.274996 520.000 18.0 
1.6500 19.636 0.250429 520.000 18.0 
1.9500 18.539 0.250003 520.000 18.0 
2.2500 17.443 0.250001 520.000 18.0 
2.5500 16.346 0.250001 520.000 18.0 
2.8500 15.250 0.250002 520.000 18.0 
3.0000 14.700 0.250002 520.000 18.0 
OIL PRODUCTION IS 3.3338 POUNDS 
WATER PRODUCTION IS 0.0153 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0076 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.9761 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 0.9872 


INSTANTANEOUS MATERIAL BALANCE FUR WATER IS 0.9052 


DISTANCE PRESSURE WATER TEMPERATURE TIME 
FEET PSTA SATURATION R HRS. 
0.0000 20.860 1-C0Q0000 660.000 18.0 
0.1500 20-812 0.518796 608.885 18.0 
0.4500 20.719 0.464119 567.152 18.0 
0.7500 20.598 0.433386 — 541.673 18.0 
1.0500 20.255 0.359360 528.789 18.0 
1.3500 19.488 0.264420 523.174 18.0 
1.6500 18.628 0.250182 521.006 18.0 
1.9500 17.762 0.250002 a Ow ae A | 18.0 
2.2500 16.890 0.250001 520.064 18.0 
2.25500 16.016 0.250002 520.012 18.0 
2.8500 15.142 0.250002 520.010 18.0 
3.20000 14.700 0.250002 520-010 18.0 
OIL PRODUCTION IS 229877 POUNDS 
WATER PRODUCTION IS 0.0137 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0465 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.9080 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 0.9886 


INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.7277 
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DISTANCE 
BEET 


228500 
3-0000 


PRESSURE 
PSIA 
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18.124 
18.032 
17.870 
17.713 
17.551 
17.361 
16.982 
16.459 
15.957 
15.456 
14.956 
14.700 


OIL PRODUCTION IS 
wATER PRODUCTION IS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 


CUMULATIVE MATERIAL BALANCE FOR WATER IS 
INSTANTANEOUS MATERTAL BALANCE FOR OIL IS 
INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 


DISTANCE 
FEET 


30000 


PRESSURE 


PSTA 


18.359 
18.353 
18.286 
18.269 
18.113 
17.889 
17.385 
16.793 
16.201 
15.603 
15.004 
14.700 


OIL PRODUCTION IS 
WATER PRODUCTION IS 


WATER 


SATURATION 


1.000000 
0.548738 
0.491520 
0.466102 
0.448451 
0.402748 
0.297437 
0.251493 
0.250016 
0.250004 
0.250003 
0.250003 


4.25586 


0.0209 


POUNDS 
POUNDS 


WATER 


SATURATION 


4.0482 
0.0186 
CUMULATIVE MATERIAL BALANCE FOR 
CUMULATIVE MATERIAL BALANCE FOR 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 
INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 


1.000000 
0.544759 
0.473977 
0.457646 
0.435225 
0.371400 
0.270977 
0.250331 
0.250006 
0.250004 
0.250003 
0.250003 


POUNDS 
POUNDS 
OIt IS 
WATER IS 


TEMPERATURE 


520.000 


1.0043 
0.8859 
1.0524 
0.4480 


TEMPERATURE 


1.0309 
0.8073 
0.9775 
0.4944 


TIME 
HRS. 


27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 


TIME 
HRS. 


27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
27.0 
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DISTANCE PRESSURE WATER TEMPERATURE TIME 


FEET _  -PSTA SATURATION R HRS. 
0.0000 17.197 1.C00000 520.000 36.0 
0.1500 17.126 0.566632 520.2000 36.0 
0.4500 16.991 0.503133 520.000 36.0 
0.7500 16.864 02474316 520-000 36.20 
1.0500 16.726 0.457186 520.000 36.0 
1.3500 16.585 0.436486 520.000 36.0 
1.6500 16.407 0.368119 520.000 36.0 
1.9500 16.046 0.268540 520.000 36.0 
2-2-2500 15.655 0.250277 520.000 36.0 
2.5500 15.273 0.250008 520.000 36.0 
2.8500 14.892 0.250005 520.000 3620 
3-0000 14.700 0.250005 520.000 3620 

OIL PRODUCTION IS 521734 POUNDS 
WATER PRODUCTION IS 0.0237 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0138 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.7607 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 1.1116 


INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.3544 


DISTANCE PRESSURE WATER TEMPERATURE TIME 
FEET PSIA SATURATION R HRS. 
0.0000 16.487 1.-C00000 660.000 36.0 
0.1500 16.2493 0.564733 624.206 36.0 
0.4500 16.439 0.473330 589.351 36.0 
0.7500 16.455 0.464989 562-738 36.0 
1.0500 16.366 02447393 544.181 36.0 
1.3500 16.312 602422002 532.719 36.0 
1.6500 16.156 0.337605 526.287 36.0 
1.9500 15.844 0.257394 522-2910 36.20 
222500 15.522 0.250084 521.2228 36.0 
22-5500 15.197 0.250007 520.405 36.0 
2.8500 14.868 0.250005 519.828 36.0 
320000 14.700 0.250005 519.2828 36.0 
OIL PRODUCTION IS 427188 POUNDS 
WATER PRODUCTION IS 0.0216 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0216 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.6990 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 0.956% 


INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.2672 
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DISTANCE PRESSURE WATER TEMPERATURE TIME 


ened PSIA SATURATION R HRS. 
0.0000 16.624 1.C00000 520.000 45.0 
0.1500 16.567 0.580891 520.000 45.0 
0.4500 16.449 0.512701 520.000 45.0 
0.7500 16.2344 0.481197 520.000 45.0 
1.0500 16.234 0.463497 520.000 45.0 
1.3500 16.110 0.447969 520.000 45.0 
1.6500 16-007 0.412629 520-000 45.0 
1.9500 15.800 0.313732 520.000 45.0 
222500 15.477 0.253042 520.000 45.0 
22-5500 15.166 0.250035 520.000 45.0 
2.8500 14.856 0.250007 520.000 4520 
3-0000 14.700 0.250007 520.000 45.0 
OIL PRODUCTION IS 5-2-6776 POUNDS 
WATER PRODUCTION IS 0.0260 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0222 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.6731 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS Lell27 


INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.2889 


DISTANCE PRESSURE WATER TEMPERATURE TIME 
FEET PSIA SATURATION R HRS. 

0.0000 15.419 1.CO00000 660.000 45.0 
0.1500 15.427 0.581523 628.169 45.0 
0.4500 15.381 0.468381 595.953 45.0 
0.7500 15.404 0.468260 570.083 45.0 
1.0500 15.349 0.448852 550. 733 45.0 
1.3500 15.357 0.432258 537.700 45.0 
1.6500 15.331 0.375970 529.666 45.0 
1.9500 15.214 0.274082 524.987 45.0 
2.2500 15.071 0.250424 5222284 45.0 
2.5500 14.927 0.250012 520.531 45.0 
2.8500 14.778 0.250007 518.624 45.0 
3-0000 14.700 0.250007 518.624 45.0 

OIL PRODUCTION IS 520600 POUNDS 

WATER PRODUCTION IS 0.0232 POUNDS 

CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0171 

CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.5967 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 0.9297 


INSTANTANEOUS MATERIAL BALANCE FOR WATER IS O.1174 
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COMPUTER OUTPUT III 


Y = NABLA = 1.500 
AXS= 7px = 0.300 feet 
kt = DT = 0.015 Hrs. 
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DISTANCE PRESSURE WATER TEMPERATURE TIME 
FEET PSIA SATURATION R HRS.» 
0.0000 252448 1.000000 520.000 9.00 
0.1500 25.066 0.482922 520.000 9.00 
0.4500 24.686 0.407252 520.000 9.00 
0.7500 23.791 0.294014 520.000 9.00 
1.0500 22.537 0.251132 520.000 9.00 
1.3500 21.330 0.250006 520.000 9.00 
1.6500 20.124 0.249998 520.000 9.00 
1.9500 18.919 0.249999 520.000 9.00 
22500 17.713 0.250000 520.000 9.00 
2.5500 16.508 0.250001 520.000 9.00 
2-8500 15.303 0.250001 520.000 9.00 
3.0000 14.700 0.250001 520.000 9.00 
OIL PRODUCTION IS 1.6944 POUNDS 
WATER PRODUCTION IS 0.0078 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL [IS 1.0146 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 1.0010 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 0.9995 
INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 1.0051 
CIS LANCE .PRESSURE WATER TEMPERATURE _ TIME 
FEET PSIA SATLRATION R HRS. 
0.COCC 24-231 1-CCCCCO 660.000 9.00 
0.-15C0 24.054 0.462056 589.177 9.00 
0.45CC 22.875 == 0.4C 2780 545.909 9.00 
0.7500 23.2244 0.291069 527.821 9.00 
1.C50C 22.163 0.250959 521.930 9.00 
1.3500 21.037 0.25CCO0l 520.391 9.00 
1-65CC 1528SC 0.249998 520.061 9.00 
1-950C 18.737 0.249999 520.001 9.00 
Seve  gilemtid _ 0.25C000 | 519.992 9.00 
2 «E500 16.43C 0.25C001 519.994 9.00 
2-85CC 15 276 0.25C002 519.998 9.00 
3-CO0C 14.70C 0.25C002 519.998 9.00 
CIL PRCCUCTION IS 1.5879 POUNDS 
WATER PROCUCTION IS  _ O0.CC73 POUNDS 
CUMLLATIVE MATERIAL BALANCE FOR OIL IS 1.0618 
CUPMLLATIVE MATERIAL BALANCE FOR WATER IS 0.9813 
INSTANTANECUS MATERIAL BALANCE FOR OIL IS 1.0206 


INSTANTANEGUS MATERIAL BALANCE FOR WATER IS 0.9813 
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DISTANCE PRESSURE WATER TEMPERATURE 
FEET PSIA SATURATION R 
0.0000 21.992 1.000000 520.000 
0.1500 21.815 0.525016 520.000 
0.4500 21.608 0.474276 520.000 
0.7500 21.299 0.437938 520.000 
1.0500 20.795 0.355678 520.000 
1.3500 19.881 0.262356 520.000 
1.6500 18.928 0.250148 520.000 
1.9500 17.989 0.250003 520.000 
2.2500 17.049 0.250002 520.000 
2.5500 16.110 0.250003 520.000 
2.8500 15.170 0.250004 520.000 
3.0000 14.700 0.250004 520.000 
OIL PRODUCTION IS 3.1353 POUNDS 
WATER PRODUCTION IS 0.0144 POUNDS 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0163 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.9256 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 1.0222 
INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.8015 
DISTANCE PRESSURE WATER TEMPERATURE 
FEET PSTA SATURATION R 
0.COCC 21.787 1.CCOCOO 660.000 
0.15CC 7p eg wat 0.530240 608.715 
0.4500 21.646 0.473011 567.159 
0.7500 21.4592 0.44C003 541.700 
1.C5CO 21.113 0.370825 528.795 
1.3500 “wa eaue, amen ~~ “Searels 
1.€5C0 19.284 0.250295 521.003 
1.95C0 18.276 0.250001 520.271 
2.25C0 17.257 0.250001 520.058 
2.5500 16.235 0.25C003 520.008 
2.85CC 15.213 0.250004 520.010 
Seca. ae OC ees 520.010 
CIL PRCCUCTICN IS 3.1891 POUNDS 
WATER PROCUCTION IS C.C146_ POUNDS etd : 
CUMULATIVE MATERIAL BALANCE FOR OIL IS 1.0351 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.9591 


-INSTANTANEGUS MATERIAL BALANCE FOR OIL [TS 2.0032 _.. 


INSTANTANECUS MATERIAL BALANCE FOR WATER IS 0.8581 
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DISTANCE PRESSURE WATER 
43m - thee LA SATURATION 
0.0000 20.792 1.c00000 
0.1500 20.664 0.548801 
024500 20.504 0.494813 
0.7500 20.292 0.469773 
1.0500 20.020 0.450660 
1.3500 19.607 0.397868 
1.6500 18.869 0.291021 
1.9500 17.915 0.251093 
222500 16.995 0.250014 
225500 16.077 0.250005 
2.8500 15.159 0.250006 
30000 14.700 0.250006 

OIL PRODUCTION IS 4.4801 POUNDS 

WATER PRODUCTION IS 0.0205 POUNDS 


CUMULATIVE MATERIAL BALANCE FOR OIL IS 
CUMULATIVE MATERIAL BALANCE FOR WATER IS 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 


520.000 


1.0195 
0.8840 
1.0258 


320.033-——__—_-# 7.00 


INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.7855 
CISTANCE PRESSURE WATER TEMPERATURE 
FEET PSIA SATURATION R 
0.<COCC 1$.16S 1.ccccao 660.000 
0.15CC 19.144 °£58247 618.368 
0.45CC 19.088 0.451312 580.697 
C275CC 1$.C24 0.463697 553.557 
Pa ot els 19.891 0.446138 536.784 
1.350C - AL8<643-— ———"6,4¢614646—-——__ _ -—_<$27-2686- 
T.650C 14.1C7 0.297491 523.248 
, 1.-S95CC 17.351 0.251512 521.248 
2eea0C 16.6C2 0.250014 520.431 
6500 15.844 0.250004 520.136 
2-85CC 15.083 0.250006 520.033 
2c ee OC: Ore ero 
CIL PRCLUCTICN IS 4.4768 POUNDS 
wATER PROCUCTICN IS C.C205 POUNDS ’ 
CUMLLATIVE MATERIAL BALANCE FOR OIL IS T.0257 
CUMLLATIVE MATERIAL BALANCE FOR WATER IS 0.8886 
INSTANTANECUS MATERIAL EALANCE FOR OTL IS SA4.C. 
INSTANTANECUS MATERIAL BALANCE FOR WATER IS 0.6425 
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DISTANCE PRESSURE WATER TEMPERATURE TIME 
FEET PSIA SATURATION R HRS. 

0.0000 19.594 1.C00000 520.000 36.00 
0.1500 19.490 0.566113 520.000 36.00 
0.4500 19.358 0.509392 520.000 36.00 
0.7500 19.182 0.481667 520.000 36.00 
1.0500 18.966 0.467165 520.000 36.00 
1.3500 18.715 0.456059 520.000 36.00 
1.6500 18.382 0.423150 520.000 36.00 
1.9500 17.802 0.328592 520.000 36.00 
2.2500 16.899 0.255221 520.000 36.00 
2.5500 16-016 0.250059 520.000 36.00 
2.8500 15.139 0.250008 520.000 36.00 
3.0000 14.700 0.250008 520.000 36.00 

OIL PRODUCTION IS 5.27692 POUNDS 

WATER PRODUCTION IS 0.0265 POUNDS 

CUMULATIVE MATERIAL BALANCE FOR OIL [5S 1.0206 

CUMULATIVE MATERIAL BALANCE FOR WATER IS 0.8544 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 1.0236 
INSTANTANEOUS MATERIAL BALANCE FOR WATER IS 0.7488 

CISTANCE PRESSURE wATER TEMPERATURE TIME 
FEET PSIA SATURATION R HRS. 
C.COCcCc 17.642 LorCCCCOO 660.000 36.00 
C.15CC 17.832 2579027 624.182 36.00 
0.4500 17.788 O.&C0609 589.796 36.00 
0.75CC 17.762 0.471177 562.888 36.00 
1.C5CC 17.7CE 0.458223 544.362 36.00 
1.25CC 17.545 ne FPF eGR TO — —_ 532.800 36-00 
1.6500 17.305 0.4C4672 526.278 36.00 
1.S5CC 16.865 0.300452 522.895 36.00 
2225CC 16.245 0.251738 521.227 36.00 
2.55CC 15.636 0.250019 520.412 36.00 
2.85CC 16.014 0.250007 519.850 36.00 
3-COCC 14.7CC = eC I BU 36-00 

CIL PRCCUCTICN IS 5.4708 POUNDS 

wATER PROCUCTICN IS C.C251 POUNDS 

CUMLLATIVE MATERIAL BALANCE FOR OIL IS 1.0219 

CUMLLATIVE MATERIAL BALANCE FOR WATER IS 0.8111 
INSTANTANEOUS MATERIAL BALANCE FOR OIL IS 1.0012 
INSTANTANECLS MATERIAL BALANCE FOR WATER [S 0.25215 
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